Equation of State, Nonlinear Elastic Response, and Anharmonic Properties of Diamond-Cubic Silicon and Germanium: First-Principles Investigation

2.1 Nonlinear Elasticity Theory

Nonlinear elasticity theory is the theory of elasticity under finite deformation. Suppose the underformed state X moves to the deformed state x when a solid body undergoes a finite deformation. We denote the undeformed coordinates of a point in configuration X by (X₁, X₂, X₃) and represent that in configuration x by (x₁, x₂, x₃). Then, the deformation gradient matrix can be defined as

(1)αij = ∂xi/∂Xj, (1)

where the subscripts i and j vary from 1 to 3. The Lagrangian strain η_ij can be expressed as [8]

(2)ηij = 12(∑kαkiαkj − δij), (2)

where δ_ij is the Kronecker delta (unity when i=j, zero otherwise). And the free energy F and internal energy U per unit mass can be expanded respectively as Taylor series in terms of the Lagrangian strain [9–11]

(3a)F(X, ηij, T) = F(X,  0, T) + V(X)∑ijTijT(X)ηij + 12!V(X) ∑ijklCijklT(X)ηijηkl + 13!V(X)∑ijklmnCijklmnT(X)ηijηklηmn + … (3a)

(3b)U(X, ηij, S) = U(X, 0, S) + V(X)∑ijTijS(X)ηij + 12!V(X) ∑ijklCijklS(X)ηijηkl + 13!V(X)∑ijklmnCijklmnS(X)ηijηklηmn + …. (3b)

Then the isothermal and adiabatic elastic constants at state X can be defined as [12]

(4a)CijklT(X) =  1V(X)∂2F∂ηij∂ηkl|Xη′mnT, CijklmnT(X) = 1V(X)∂2F∂ηij∂ηkl∂ηmn|Xη′pqT, (4a)

and

(4b)CijklS(X) =  1V(X)∂2U∂ηij∂ηkl|Xη′mnS, CijklmnS(X) = 1V(X)∂2U∂ηij∂ηkl∂ηmn|Xη′pqS, (4b)

where V(X) is the volume of the system at state X, and the subscripts η′mn denotes that the elements of the strain tensor other than those involved in the partial derivative are held constant. T and S represent the isothermal and adiabatic processes, respectively. T_ij(X) appeared in formulas (3a) and (3b) is the component of the stress tensor applied on the state X, which can be defined as

(5)TijT(X) = 1V(X)∂F∂ηij|X η′mn T, TijS(X) = 1V(X)∂U∂ηij|X η′mn S. (5)

Specifically, they can be expressed as

(6a)TijS = TijT = −Pδij hydrostaticpressure, (6a)

(6b)TijS = TijT = −Ptitj uniaxialstressalong t, (6b)

where t is a unit vector. Here, we do not distinguish the two types of elastic constants, since the first-principles calculations are performed at 0 K, i.e., F=U − TS=U and C^S=C^T. So is T_ij(X).

If we introduce the density ρ₀ of the unstrained crystal and apply the Voigt convention (η₁₁→η₁, η₂₂→η₂, η₃₃→η₃, η₂₃→η₄, η₁₃→η₅, and η₁₂→η₆), (3a) can be simplified as

(7)ρ0F(X,ηij) = ρ0F(X, 0) + ∑iTi(X)ηi + 12!∑ijCij(X)ηiηj  + 13!∑ijkCijk(X)ηiηjηk + ⋯. (7)

For cubic crystals, there are three independent second-order elastic constants C₁₁, C₁₂, and C₄₄ and six independent third-order elastic constants C₁₁₁, C₁₁₂, C₁₂₃, C₁₄₄, C₁₆₆, and C₄₅₆. Thus, (7) can be further written as

(8)ρ0[F(X, ηij) − ρ0F(X, 0)] = ∑iTi(X)ηi + φ2 + φ3, (8)

where

(9)φ2 = 12C11(η12 + η22 + η32) + C12(η1η2 + η2η3 + η1η3) + 12C44(η42 + η52 + η62) (9)

and

(10)φ3 = 16C111(η13 + η23 + η33)+12C112(η12η2 + η12η3 + η22η1 + η22η3 + η32η1 + η32η2) + C123η1η2η3 + 12C144(η42η1 + η52η2 + η62η3) + 12C166(η42η2 + η42η3 + η52η1 + η52η3 + η62η1 + η62η2) + C456η4η5η6. (10)

Using (8)–(10), we can obtain the elastic constants by applying a series of specific homogeneous deformations on a crystal and calculating the total energies of the undeformed and deformed states, respectively. In the following, the special deformation types and calculation details of the total energies are presented.

2.2 Deformation Types

Since the single Si and Ge at 0 K and 0 GPa belong to the cubic structure, and there are six independent third-order elastic constants, six special deformation types at the least are needed to calculate their complete set of the third-order elastic constants. Moreover, the six deformation types should be as simple as possible to reduce the number of the strain components that will appear in (9) and (10). Following this principle, we select the six deformation types [3] η_A, η_B, η_C, η_D, η_E, and η_F shown in Table 1. Then, we obtain the corresponding energy expressions f_α by inserting the six different deformation strains into (8)–(10),

(11)fα = ρ0[F(X, η) − ρ0F(X, 0)] = P2η2 + P3η3 + o(η4), (11)

where f_α=f_A, f_B, f_C, f_D, f_E, and f_F; the coefficients P₂ and P₃ are related to the second- and third-order elastic constants; in practice, the strain magnitude η is set from −|η_max| to +|η_max| in steps of 0.0005. Here, it should be pointed out that the term ∑ijTij(X)ηij does not appear in (11), for T_ij(X) representing that the stress applied on the reference state is zero at 0 GPa.

For each deformation type listed in Table 1, the deformed state can be obtained by using the equation of r′i = αijrj, where r′i and r_j represent respectively the deformed and undeformed lattice vectors; α_ij can be obtained by the following formula:

(12)αij = δij + ηij − 12∑kηkiηkj + 12∑rkηrkηriηkj − 58∑kmnηkjηmkηmnηni + ⋯. (12)

In general, for a given η, α_ij is not unique, but this is not a problem since the Lagrange strain brings the rotational invariance of total energy [9]. And for a system without rotation, (12) will provide a unique relation.

2.3 Computational Details

2.3.2 Determination of the Third-Order Elastic Constants

Since the third-order elastic constants are very sensitive to the k-points, the maximum strain magnitude, and the order of the polynomial used to fit (11), we take Si as an example to test the dependence of the third-order elastic constants on the three parameters. First of all, we investigate the dependence of the third-order elastic constants on the k-points and display it in Figure 1. As shown in Figure 1, all the third-order elastic constants of Si converge after the k-point mesh size reaches 15×15×15. And the relative difference between two successive values of an examined constant is <1.0 %, which is exactitude enough for the calculations of the third-order elastic constants. Such a test is also performed for Ge, and results show that EcutoffGe = 650 eV and 17×17×17 are sufficient to calculate the third-order elastic constants. Subsequently, the dependence of the third-order elastic constants on the maximum strain magnitude |η_max| is tested and illustrated in Figure 2. From Figure 2, it can be clearly found that the third-order elastic constants converge well after the maximum strain parameters |η_max| up to 0.009. Therefore, |η_max|=0.009 is employed to calculate the third-order elastic constants of Si. Similarly, |η_max|=0.015 is determined for Ge. Based on the parameters determined above, we obtain the relations of η and f_α(η) (α=A, B, C, D, E, and F). Then the remaining step is to fit these relationships using proper polynomials. For ascertaining the suitable order of the polynomials, the residuals obtained respectively from the third- and fourth-order polynomial fits are compared in Figure 3a–f. Figure 3a–f show that the residuals obtained from the fourth-order polynomial fit are clearly smaller than those obtained from the third-order polynomial fit. Thus, the fourth-order polynomial is employed to fit the relations of f_α(η) and η. Then, we present the f_α(η) ∼ η data of Si in Figure 4 and fit them using the fourth-order polynomial. From Figure 4, the second- and third-order elastic constants of Si are obtained. In the same way, the elastic constants of the Ge are obtained.

Figure 1:

Dependence of the third-order elastic constants C₁₁₁ and C₁₁₂ (a), C₁₂₃ and C₁₄₄ (b), and C₁₆₆ and C₄₅₆ (c) on the k-points (energy cutoff of 650 eV is applied for all points).

Figure 2:

Convergence tests of the third-order elastic constants of Si as a function of the maximum strain |η_max|.

Figure 3:

Comparisons of the residuals obtained respectively from the third- and fourth-order polynomial fit. (a), (b), (c), (d), (e), and (f) represent respectively the relations between the residual of f_A, f_B, f_C, f_D, f_E, f_F and η.

Figure 4:

The strain energy density as a function of the Lagrangian strain of Si, where the discrete points and the solid lines represent respectively the results obtained from the first-principles calculations and the nonlinear theory; f_A, f_B, f_C, f_D, f_E, and f_F denote different energy expressions corresponding to various deformation modes.

3.1 Equation of State

Within the framework of the LDA, we calculate a series of total energies of Si and Ge in both diamond-cubic and β-tin structures and fit them using the fourth-order Brich–Murnaghan equation [33]

(13)E(V) = ∑n = 14anV−2n/3. (13)

Then, the equation of states of Si and Ge are obtained and displayed in Figure 5. It is known that the diamond-cubic structure will coexist with the β-tin phase at the pressure where the Gibbs free energies, G=E + PV−TS, of the two structures are equal. At zero temperature, the Gibbs free energy is equal to the enthalpy H=E + PV, where the pressure is given by the identity P=−∂E/∂V. That is, the phase transition pressures from the diamond-cubic structure to the β-tin phase can be determined by using the enthalpy method. Therefore, we display their enthalpy difference versus pressure in Figure 5.

Figure 5:

Energy versus primitive volume for Si and Ge in the diamond-cubic and β-tin structures. Inset: enthalpy differences ΔH of the β-tin structures relative to diamond phases as a function of the pressure.

Figure 5 shows that the phase transition (diamond cubic → β-tin) of Si occurs at about 7.2 GPa, which is exactly between the previous theoretical data 7.0 [34] and 8.4 GPa [20, 35, 36], but less than the available experimental values 11.7 GPa [37] and 10.0–12.5 GPa [38, 39]. For Ge, the corresponding phase transformation is near 9.7 GPa, which is in general agreement with the previous theoretical results, 9 GPa [40], and also lower than the existing experimental results, 10.5 ± 0.2 [41] and 8.1 ± 0.3–10.6 ± 0.5 GPa [29]. Here, it is worth pointing out that such underestimation of the phase transition pressure is typical in the LDA calculations.

3.2 Results of the Third-Order Elastic Constants

The second- and third-order elastic constants of the diamond-cubic Si and Ge are listed in Table 3 and compared with other theoretical data [3, 4, 19, 42] and the available experimental results [27, 43–48]. For the second-order elastic constants, it is obvious that the present data for both Si and Ge are well consistent with the existing experimental results [27, 43, 45, 46] and foregoing theoretical values [3, 4, 19]. The differences do not exceed the typical error (10 %), which is considered the typical error in elasticity calculations based on the DFT. Such agreement suggests that the elasticity calculations are performed at high accuracy and the elastic constants obtained in the present work are reliable. Here, it should be pointed out that some elastic constants can be determined by several fits, while the obtained values are slightly different (e.g., for Si, C₄₄=77.06, 77.33, and 78.20 from f_D, f_E, and f_F, respectively). In such cases, the average of these data is given in Table 3. In addition, Table 3 clearly shows that the second-order elastic constants are in better agreement with the available experimental results [27, 43, 45, 46] than those of other theoretical values [3, 4, 19].

As regards the third-order elastic constants of Si and Ge, overall agreement of our calculated data and available experimental [27, 43–48] and other theoretical results [3, 4, 19, 42] can be seen from Table 3. The discrepancies between our calculated values and experimental results may be derived from the following two aspects. First, the experimental measurements of the third-order elastic constants are still difficult to date; i.e., the reported experimental results [27, 43–48] are determined with considerable uncertainties and even exhibit a difference between different groups (e.g., [43] and [44] shown in Table 3). Second, theoretical results of the third-order elastic constants are obtained at 0 K, while the experimental values are often determined in conditions which are far from the ideal case of 0 K. Particularly, for the cases of Si and Ge, the temperature effects on the third-order elastic constants can be identified and should not be ignored (see [31, 44, 49]). Moreover, the distinctions between the values obtained in the present work and other theoretical data, in our opinion, are understandable, since the theoretical results obtained at different times by different numerical procedures often have different and unknown accuracies.

4.1 Ascertainment of the Pressure Derivatives of the Second-Order Elastic Constants

For describing the nonlinear elastic properties of a material under large hydrostatic pressure, we need to introduce the effective elastic constants C_ij(P), which can be obtained by the following formula [10]:

(14)Cij(P) ≈ Cij + C′ijP + ⋯, (14)

where C′ij denotes the pressure derivative of C_ij. If C′ij is known, the effective elastic constants C_ij(P) of a material under high pressure can be calculated readily through (14). For cubic crystals, the pressure derivatives C′ij can be expressed as [10]

(15)C′11 = −2C112 + C111 + 2C12 + 2C112C12 + C11C′12 = −C123 + 2C112 − C12 − C112C12 + C11C′44 = −2C166 + C144 + C44 + 2C12 + C112C12 + C11. (15)

Using (15) and the results of the second- and third-order elastic constants, we calculate the pressure derivatives C′ij of Si and Ge. To compare them with other theoretical data [3, 50] and experimental results [49, 51], we list them in Table 4.

With regard to the diamond-cubic Si, Table 4 shows that our calculated pressure derivatives are in reasonable agreement with the existing experimental results [49, 51]. The discrepancies between our data and experimental values respectively are 9.2 %–12.9 % for C′11, 1.0 %–3.3 % for C′12, and 20 %–25 % for C′44. Though the difference between the present result of C′44 and the experimental data is large, it is much smaller than the discrepancies (64 %–66 % for [3] and 37.5 %–46.7 % for [50]) between the experimental data and the values offered in [3, 50]. We guess that the difference between theoretical data and experimental values may be originated from the temperature effects which have a large influence on the elastic constants of Si. For the pressure derivatives of the diamond-cubic Ge, however, there are no experimental and other theoretical values to date unfortunately. So it is hoped that the results obtained in this work would serve as a reference for future studies.

4.2 Calculations of the Grüneisen Constants of Long-Wavelength Acoustic Modes

The Grüneisen constant is a most important parameter for describing the anharmonic behaviour of solids such as thermal expansivity, shock deformation, and ultrasonic attenuation. In the quasiharmonic approximation, the generalized Grüneisen constants can be defined as [52]

(16)γ(q, j) = −∂lnω(q, j)∂lnV = −Vω(q, j)∂ω(q, j)∂V, (16)

where V is the crystal volume, q denotes the propagation direction, j represents the polarization vector, and ω stands for the frequency of the ith mode. In addition, Brugger defined the strain Grüneisen tensor as [53]

(17)γαβ(q, j) = −1ω(q, j)∂ω(q, j)∂ηαβ (α, β = x, y, z), (17)

where η_αβ is the Lagrangian strain. For a cubic crystal, the generalized Grüneisen constant is equal to one third of the trace of the strain Grüneisen tensor. Namely,

(18)γ(q, j) = 13[γxx(q, j) + γyy(q, j) + γzz(q, j)]. (18)

According to the motion equation [9] of a long-wavelength elastic wave, the element of the strain Grüneisen tensor for a given q can be written as

(19)γαβ(q, jj′) = −14ρ0vqjvqj′∑μνςξ(Sμν|ςξ|αβ + Sμν|ςξ|βα)  × qνqξeμ(qj)eς(qj′), (19)

where ρ₀ is the density of the undeformed crystal, v represents the sound velocity, e(q, j) denotes the polarization vector of branch j and propagation direction q, and S is the combination of the second- and third-order elastic constants

(20)Sαμ|βν|ςξ = Cαμ|βν|ςξ + δαβCμν|ςξ + δαςCβν|μξ + δβςCαμ|νξ. (20)

For a specific q and j, γ(q, j) can be expressed in terms of the second- and third-order elastic constants (see the appendix of [52]) by combining (18)–(20) and a series of equations listed in Tables 1–3 of [8]. Using these equations displayed in the appendix of [52] and our results of the elastic constants, we compute the thermodynamic Grüneisen constants of Si and Ge and list them in Table 5, together with other theoretical data [3, 52] and available experimental results [44, 46]. With regard to Si, the comparisons show that our calculated results are in better agreement with the experimental values [44, 46] than those of other theoretical data [3, 52]. This indicates that our results of the elastic constants are reliable and can be used to further calculate the ultrasonic nonlinear parameters. In addition, this phenomenon also illuminates that the Grüneisen constants are very sensitive to the values of the elastic constants, since our elastic constants are also in better agreement with the experimental data than that of the other theoretical calculations. To the best of our knowledge, there are no experimental and theoretical data on the Grüneisen constants of Ge till now. Hence, the Grüneisen constants of Ge calculated in this work may be severed as a reference for future studies.

4.3 Determination of the Ultrasonic Nonlinear Parameters

The ultrasonic nonlinear parameters, obtained by the formula [9]

(21)β(q,jj′) = −1ρvqj2∑αβμνςξSαμ|βν|ςξqμqνqξ × eα(q, j′)eβ(q, j)eς(q, j), (21)

are the indispensable quantities in the second-harmonic generation experiments. Combining with (19), we obtain the general relation between the ultrasonic nonlinear parameters and the Grüneisen constants, which is written as

(22)β(q, jj′) = 2∑αβγαβ(q,jj′)qαeβ(qj). (22)

For the longitudinal (L) modes along the main-symmetry direction of cubic crystals, (22) can be reduced as the following three formulas:

(23)β([1, 0, 0], LL) = 2γxx([ε, 0, 0], L),β(12[1, 1, 0], LL) = 2{γxx([ε, ε, 0], L) + γxy([ε, ε, 0], L)}andβ(13[1, 1, 1], LL) = 2{γxx([ε, ε, ε], L) + 2γxy([ε, ε, ε],L)}. (23)

Thus, the three ultrasonic nonlinear parameters can be explicitly expressed in terms of the second- and third-order elastic constants by using the equations given in the appendix of [52]. Then we calculate the three ultrasonic nonlinear parameters of Si and Ge by using the results of the second- and third-order elastic constants and list them in Table 5. As displayed in Table 5, three ultrasonic nonlinear parameters of Si are in reasonable agreement with the existing experimental results [44] measured at 300 K, and the discrepancies for β₁, β₂, and β₃, respectively, are 7.9 %, 2.8 %, and 15.1 %. It is a pity that there are no experimental and theoretical results on the ultrasonic nonlinear parameters of the Ge up to now. Actually, the comparison between our results and experimental data as well as other theoretical values is not our aim. Instead, our aim is to illuminate some applications of the third-order elastic constants and introduce the way of calculating the ultrasonic nonlinear parameters on the basis of the elastic constants.