Density Functional Study of the Carbon Dependence of the Structural, Mechanic, Thermodynamic, and Dynamic Properties of SiC Alloys

Using a density functional scheme, for the first time the carbon dependence on the structural, dynamic, thermodynamic, and dynamic properties of $$\hbox {Si}_{1-x}\hbox {C}_{x} $$Si1-xCx alloys ($$x=0.0$$x=0.0 to 1.0 in steps of 0.125) has been investigated. The structural properties of these materials, in particular, the composition dependence of the lattice parameter and bulk modulus, are in excellent agreement with experimental data and follow a quadratic law in (x). A nonlinear relationship is found between the elastic constants $$C_{11}$$C11, $$C_{12}$$C12, and $$C_{44}$$C44 and the carbon concentration (x). The behavior of the acoustical and optical phonon frequencies at high-symmetry points $$\varGamma $$Γ, X, and L is predicted. Through the quasi-harmonic Debye model, in which the photonic effects are taken into account, the Debye temperature, the heat capacity, the Helmholtz free energy, the internal energy, and the entropy are determined for the $$\hbox {Si}_{1-x}\hbox {C}_{x }$$Si1-xCx compounds.


Introduction
In view of their technological promise, Si 1−x C x alloys have attracted increasing research interest in the last few years. Although they are members of the binary group IV alloy family, they have remarkable properties, which set them apart from other binary alloys, such as Si 1−x Ge x [1,2]. These alloys containing carbon exhibit both B H. Langueur langueur_hassina@yahoo.fr 1 Optoelectronic and Compounds Laboratory, Department of Physics, FERHAT ABBAS University, 19000 Sétif, Algeria physical and chemical properties that differ significantly from the general trends of the group IV family.
Silicon and carbon, both crystallizing in the zincblende structure, form a continuous series of Si 1−x C x solids with composition x ranging from 0 to 1, where x denotes the mole fraction of C. By changing the alloy composition, the physical properties can be controlled at arbitrary values of composition x between those of Si and C, and in some cases, they can be quite different from those of the constituent materials. To this purpose, the effect of composition on physical properties of Si 1−x C x assumes particular importance.
In this work, we are interested in the prediction of the effect of the carbon composition on the structural, elastic, thermodynamic, and dynamic properties of the zincblende Si 1−x C x compounds, and we present for the first time a complete and comparative study of these properties for several carbon concentrations from 0 to 1 (in steps of 0.125).
The paper is organized as follows: Sect. 2 contains a brief description of the method of calculation. Section 3 contains our results compared to the available experimental and theoretical data. Section 4 contains the conclusions.

Method of Calculation
We have performed first-principle total-energy calculations using the pseudopotential method based on the density functional theory (DFT) [3,4], within the local-density approximation. We employed norm-conserving Troullier-Martins pseudopotentials [5]. The exchange-correlation energy of electrons is described in the local-density approximation with the Teter-Pade parameterization [6]. The calculations were carried out using the ABINIT code [7].
The Kohn-Sham (KS) single-particle functions were expanded in a plane-wave basis set. Self-consistent solutions of the (KS) equations were obtained by sampling the irreducible Brillouin zone with the special k-point method. Well-converged results were obtained using a kinetic energy cutoff of 70 Hartree and with a set of 28 kspecial points, which correspond to 6 × 6 × 6 k-points mesh in the Monkhorst-Pack notation [8]. Having obtained self-consistent solutions of the (KS) equations, phonon frequencies and elastic constants are obtained using the self-density functional perturbation theory (DFPT) [9] which allows the calculations of the dynamic matrix at arbitrary q vectors. We ensure the convergence of the phonon frequencies to 2 cm −1 to 3 cm −1 .
To model the alloy at different compositions (x), we have used the virtual crystal approximation (VCA) [10].VCA treats an alloy as a perfect periodic crystal, whose properties are modeled as a linear average of those of Si and C atoms. Elemental ionic pseudopotentials of Si and C are combined to construct the virtual pseudopotential of the Si 1−x C x :

Structural Properties
In the present work, the structural properties of the binary compounds Si and C are calculated in the zincblende structure. Then, the semiconductor binary alloy is of the type Si 1−x C x . We have started our pseudopotential calculation of the structural properties with the zincblende structure, and let the calculation forces move the atoms to their equilibrium positions. For the considered structures, we performed the structural optimization by calculating the total energies for different volumes around the equilibrium cell volume V 0 of the Si and C compounds and their alloys. The calculated total energies are fitted to Murnaghan's equation of state [11] to determine the ground-state properties as the equilibrium lattice constant a 0 , the bulk modulus B 0 , and its pressure derivative B . The equilibrium lattice parameters have been computed for several concentrations (x) by minimizing the crystal total energy calculated for different values of the lattice constant. The initial lattice parameters were taken from previous experimental and theoretical values as a starting point for geometry optimization. The calculated lattice parameters are listed in Table 1.
The calculated lattice constants for Si and C agree quite well with those experimentally obtained in [12], and they differ by less than 0.77 % and 0.75 %, respectively. It is known that the LDA calculation generally gives an ∼1 % underestimation of the lattice constants. The theoretical values by Bernard and Zunger [13] and Sze [16] are similar for Si and C (respectively) with ours.
In Fig. 1, we present our VCA (in steps of 0.125) calculated lattice constant as a function of C concentration, which is fitted by the following quadratic equation: One can see from Table 1 that like for Si and C, the LDA underestimates the lattice parameter of Si 1−x C x . Our calculated bulk modulus agrees well with the experimental values [17,18] for Si and C, respectively, and as a function of the carbon concentration, it is determined by the following equation:

Elastic Constants
It is very important to investigate the elastic properties of materials because they provide a link between the mechanical and dynamic behavior of crystals and give important information concerning the nature of the forces operating in solids. In particular, they provide information on the stability and stiffness of materials. Values of elastic constants give significant information about the bonding characteristics between adjacent atomic planes and the anisotropic character of the bonding and struc-  tural stability [22]. Elastic properties of a cubic single crystal are described by three independent constants C 11 , C 12 , and C 44 . The combination of these three constants gives more information and data, such as acoustic wave velocities, elastic moduli for crystals in polycrystalline states (bulk modulus B, shear modulus G, Young's modulus E, Poisson's ratio σ ), and some thermodynamic properties such as the Debye temperature.
Here, we have computed the elastic response of the Si 1−x C x alloys with respect to the strain perturbations, obtaining the second derivatives of the total energy with all the perturbations in order to obtain the elastic constants as implemented in the ABINIT code.
The obtained single-crystal elastic constants C 11 , C 12 , and C 44 for Si, C, and their binary alloys in the zincblende phase are listed in Table 1. There is no experimental study in the literature about the effect of carbon concentration on the elastic constants of the binary alloys Si 1−x C x ; hence, there is no possibility for comparison with experiments, except for Si and C which are in reasonable agreement with previously reported results.
It is known that the mechanical stability conditions in cubic crystals on the elastic constants are given by C 11 − C 12 > 0, C 11 > 0, C 44 > 0, C 11 + 2C 12 > 0, and C 12 < B < C 11 . Our elastic constants in Table 1 satisfy these stability conditions, and hence, we can say that Si 1−x C x alloys in the zincblende structure are elastically stable for all carbon concentrations.
The dependence of the elastic constants C 11 , C 12 , and C 44 on the carbon concentration (x) in the VCA (in steps of 0.125) is also illustrated in Fig. 2 and described by the following expressions [in GPa]: We can see from Fig. 2 that the three elastic constants C 11 , C 12 , and C 44 increase with increasing carbon concentration x.
Once the elastic moduli are calculated, we derived other important mechanical properties (see Table 1), such as the Zener anisotropy factor A z , the isotropic shear modulus G, the Poisson ratio σ , and the Young's modulus E.
The Zener ratio is defined by [23] A z = 2C 44 / (C 11 − C 12 ) For an isotropic crystal C 11 −C 12 = 2C 44 , so A Z = 1. The magnitude of the deviation from 1 is a measure of the degree of elastic anisotropy possessed by the crystal. If A Z < 1, the crystal is stiffest along the 100 cube axes, and when A Z > 1, it is stiffest along the 111 body diagonals [24]. From the computed A Z values in Table 1, A Z is found to be lower than 1 for all C concentrations of Si 1−x C x alloys, which means that they are stiffest along the 100 body diagonals. We can also see from Table 1 that the Si 1−x C x alloys became more anisotropic with decreasing concentration x of C.
Our results for the Zener ratio are fitted by the following quadratic equation: We have calculated Poisson's ratio (σ ) and Young's modulus (E) in terms of the computed data given in [25], in which the isotropic shear modulus (G) is calculated using the Voigt-Reuss-Hill homogenization method given in [26].
The value of Poisson's ratio which quantifies the stability of the crystal against shear is small (σ = 0.1) for covalent materials, whereas for ionic materials it is usually close to 0.25 [27].
Our calculated Poisson's ratios decrease from 0.22 to 0.08 with increasing concentration x. Therefore, the covalent contribution in inter-atomic bonding became higher with increasing C concentration x, and we can say that these materials became mechanically more stable with decreasing carbon concentration.
The Young's modulus (E), which is defined as the ratio of stress to strain and is used to provide a measure of the stiffness of the solid, is found to increase with increasing concentration, so we can say that these materials became much harder and more stiffer going from x = 0 to x = 1.
Our calculated Poisson's ratio (σ ), Young's modulus (E), and the isotropic shear modulus (G) versus concentration (x) in the VCA (in steps of 0.125) are described by the following expressions: There is no study in the literature about the effect of the carbon concentration on these parameters, so our calculated results can be considered as predictions.

Debye Temperature
One of the most important parameters that determines the thermal characteristics of materials is the Debye temperature θ D . As a rule of thumb, a higher θ D implies a higher associated thermal conductivity and melting temperature. The knowledge of such a numerical figure is essential for developing and manufacturing electronic devices [24,27]. The Debye temperature (θ D ) of the considered compounds is estimated from the average sound velocity ν m in terms of the following equation [28]: where h is the Planck's constant, k B is the Boltzmann's constant, n is the number of atom per molecule, N A is the Avogadro's number, ρ is the density, and M is the molecular weight. ν m is the average wave velocity and it is dependent on the transverse (ν t ) and the longitudinal (ν l ) wave velocities according to the following formulas [28]: ν l and ν t can be estimated from the shear modulus G and the bulk modulus B by using Navier's equation as follows [28]: The calculated values of average, transverse, and longitudinal wave velocities and the Debye temperature for Si 1−x C x alloys are presented in Table 2.
For the Debye temperature, we have found that θ D increases with increasing carbon concentration (x). The increase in the Debye temperature can be attributed to the progressive increase in the average sound velocity ν m .
Our result for the carbon concentration effect on the Debye temperature is given in Table 2 and is correlated by the following equation: θ D = 633.42 + 208.88x + 975.06667x 2 − 1980.16x 3 + 2408.53333x 4 (15) To our knowledge, there are no reported results of the effect of carbon concentration on average, transverse, and longitudinal wave velocities of the Si 1−x C x alloys; so there is no possibility for comparison.

Temperature Dependence of the Thermodynamic Quantities
Using the concept of the phonon, we calculate the thermal properties of the crystal, in particular, the phonon contribution to the Helmholtz free energy (F), the internal energy (E), the entropy (S), and the heat capacity (C). We calculate the temperature dependence of these parameters at a constant volume of Si 1−x C x alloys for various concentrations (x) in the VCA. The results are plotted in Fig. 3. Figure 3a displays the internal energy of Si 1−x C x alloys versus (x) as a function of the temperature. Our results suggest that, above 300 K, the total internal energy increases almost linearly with temperature for all concentrations (x = 0, 0.25, 0.500.75, 1.0). At high temperatures, the internal energy tends to display k B T behavior. Figure 3b shows the variation of the free energy versus (x) as a function of the temperature. Overall profiles of all plots show similar characteristics, and the free energy decreases gradually with increasing temperature.
The zero-temperature values F 0 and E 0 can be calculated from the following expression, generated by the zero-point motion [29]: where n is the number of atoms per unit cell, N is the number of unit cells, ω max is the largest phonon frequency, ω is the phonon frequency, and g(ω) is the normalized phonon density of states with The calculated F 0 = E 0 = (13.021, 14.38, 16.69, 21.45 and 33.85) kJ·mol −1 for x = 0, 0.25, 0.50, 0.75 and 1, respectively. It can be clearly seen that the free energy and the internal energy increase with increasing carbon concentration (x). The variation of the entropy with temperature is given in Fig. 3c for the same temperature range and for the same concentrations (x = 0, 0.25, 0.500.75, 1.0). Here the entropy increases rapidly with increasing temperature and decreases with increasing carbon concentration (x).
In Fig. 3d, we show the temperature dependence of the total heat capacity (C V ). One can see that C V shows the same behavior in the temperature range from 0 to 50 K, for all concentrations (x= 0, 0.25, 0.50 0.75, 1.0). For T < 600 K, C v increases rapidly with the temperature and with decreasing carbon concentration, and for T > 600 K, C v increases slowly with the temperature and approaches a constant called the Dulong-Petit limit.
The effect of the carbon content on the Helmholtz free energy (F), the internal energy (E), the entropy (S), and the heat capacity (C) is fitted by the following equations:  Fig. 3 Temperature dependence of internal energy, free energy, entropy, and heat capacity for Si 1−x C x compounds F kJ·mol −1 = 6.65 + 11.28x − 20.87x 2 + 37.84x 3 (18) E kJ·mol −1 = 18.39 + 9.52x − 32.88x 2 + 41.28x 3 (19) S J·mol −1 ·K −1 = 39.14 − 5.88x − 40.02x 2 + 11.45x 3 (20) We note that there are no results reported on the internal energy, free energy, and the entropy; except for the heat capacity for Si and C at T = 300 K, which are equal to 40. 26 and 12.36 H [35], respectively, and agree well with our results (39.72 and 12.33).

Phonons Dispersion
The calculation of vibrational properties is performed using DFPT in the plane-wave pseudopotential method as implemented in the ABINIT code, with the above-described  Table 3 Phonon frequencies calculated at the high-symmetry points Γ , X , and L for Si 1−x C x compounds considered in this work (cm −1 ), the related experimental data, and the other theoretical works virtual crystal approximation. Phonon band structures of Si 1−x C x have been calculated for 0 < x < 1 in steps of 0.125 of the alloy (x). Our results for the bulk phonon dispersions along several symmetry lines, for several concentrations (x) of C (x = 0, 0.25, 0.50, 0.75, and 1), are displayed in Fig. 4. When comparing phonon dispersion curves of Si 1−x C x (0 < x < 1), it becomes apparent that the general shapes in the dispersion curves are similar among these materials, and there is no splitting at the Γ point between the longitudinal and transversal optical phonon frequencies.
The optical and acoustic regions have crossed for Si, C, and Si 1−x C x alloys at the X symmetry point. The separation between the longitudinal and transversal modes for the acoustic region at the X point increases across Si → C; however, the optical region decreases with increasing carbon content.
Some numerical values at the high-symmetry points Γ , X , and L are also reported and compared with experimental and theoretical results in Table 3. Our results for Si and C agree well with experiment [36] and other theoretical [37,38] results. Note that we have no experimental data to compare with our results for Si 1−x C x alloys. Our calculations are a reliable prediction of the phonon dispersions. The composition dependence of the phonon frequencies calculated at the high-symmetry points Γ , X , and L for Si 1−x C x alloys is plotted in Figs

Conclusion
We have performed first-principle calculations in the VCA method of the structural, elastic, thermodynamic, and vibrational properties of the bulk Si 1−x C x alloy, within the LDA to the DFT using the pseudopotential plane-wave method. The calculated lattice constants of the bulk Si 1−x C x alloy are described by a quadratic function.
The elastic constants, the Zener ratio, the Young's modulus, the Poisson's ratio, the Debye temperature, and the heat capacity gave a nonlinear relationship with the C composition x.
The general shapes in the dispersion curves are similar between these materials, and there is no splitting at the Γ point between the longitudinal and transverse phonon frequencies.
The variation of the calculated optical and acoustic phonon frequencies with composition x has a nonlinear form.
For the two limiting cases, Si and C, in most cases calculations are in good agreement with available theoretical and experimental data.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.