First-Principles Study of Electronic and Elastic Properties of Hexagonal Layered Crystal MoS₂ Under Pressure

3.1 Crystal Structure and Equation of States

At ambient conditions, MoS₂ is a hexagonal layered crystal with the space group of P63/mmc (No. 194). Mo atom is located at the 2c Wyckoff position (0.3333, 0.6667, 0.250), and S atom is located at the 4f Wyckoff position (0.3333, 0.6667, 0.623) [23]. The hexagonal layered crystal structure of MoS₂ is illustrated in Figure 1, which shows that there are two MoS₂ molecules with the dihexagonal dipyramidal shape in the unit cell. In order to determine the ground state of MoS₂, we adopt the following procedure: Firstly, we select a series of different axial ratios c/a near the experimental ratio of hexagonal structure. For a given ratio c/a, a series of different lattice parameters a and c are set to calculate the total energies E and the corresponding primitive cell volumes V, and then obtained the lowest energy E_min for the given ratio c/a. This procedure is repeated over a wide range of c/a. Finally, by fitting the E_min−V data to the Vinet Murnaghan equation of state (EOS) [24],

Figure 1:

The hexagonal structure of MoS₂: (a) side view and (b) top view, where green and yellow spheres in the sketch represent the Mo atoms and S atoms, respectively. The grey, green, and magenta arrows represent the X, Y, and Z directions, respectively.

(1)ln(Px23(1 − x)) = lnB0 + a(1 − x), x = (VV0)1/3, (1)

where V=V(0, T) is the zero-pressure equilibrium volume, derived by integration of the thermodynamic definition of the thermal expansion coefficient α(T)=V⁻¹(∂V/∂T)

(2)V(0, T) = V(0, 0)exp∫0Tα(T)dT. (2)

B₀ and a( = 3(B′0 − 1)/2) are the fitting parameters, and we can obtain the equilibrium structural parameters a₀, c₀, c₀/a₀ and equilibrium primitive cell volume V₀, bulk modulus B₀ of MoS₂ at P=0 GPa and T=0 K. Our lattice parameters (a₀=3.1393 Å, c₀=12.4318 Å, c₀/a₀=3.96) of MoS₂ are consistent with the experimental values (a_exp=3.169 Å, c_exp=12.324 Å, c_exp/a_exp=3.8889) [23] with the errors of −0.94 %, 0.88 %, and 1.8 %, respectively. The obtained equilibrium cell volume (V₀=106.1 ų) of MoS₂ deviates from the experimental data (V=107.18 ų) by about −1.0 %, which is consistent with the fitting result of Chi et al. (V=106.38 ± 0.30 ų) [25]. It ensures the reliability of our calculations. The obtained bulk modulus B₀ of MoS₂ at 0 GPa and 0 K is 50.86 GPa, which is a little lower than the calculated value of 63.36 GPa by Li et al. [26], but it is consistent well with the X-ray diffraction (XRD) experiment value of 47.65 ± 0.30 GPa [25].

In Figure 2, we illustrate the pressure dependences of the normalised lattice parameters a/a₀, c/c₀, and the normalised cell volume V/V₀ of MoS₂, where a₀, c₀, and V₀ are the zero pressure equilibrium lattice parameters, together with the experimental results by Nayak et al. in their XRD experiment using a symmetric DAC [27]. Our results show that below 5 GPa the rate of reduction of lattice parameters with pressure is different along in-plane and cross-plane directions: the equilibrium ratio c/c₀ decreases more quickly than a/a₀, which indicates that the compressibility along c axis is larger than along a axis. This anisotropy vanishes above 25 GPa, which is well consistent with the experimental data [28] and the theoretic result by Guo et al. from the projector augmented wave method [29]. Our calculated lattice parameters are also in agreement with the experimental data of Bandaru et al. in their X-ray diffraction experiments [30]. Also, in Figure 2, it can be seen that as the pressure increases from 0 to 25 GPa, the c/a ratio decreases by 8.08 %. The basal plane linear compressibility, d(lnc)/dP=−0.007 GPa⁻¹ is about 3.5 times as large as the interlayer value d(lna)/dP=−0.002 GPa⁻¹, which suggests that lattice constant c is more sensitive to pressure than lattice constant a.

Figure 2:

The normalised parameters a/a₀, c/c₀, V/V₀, and ratio c/a of MoS₂ at T=0 K as a function of pressure. Our results are presented by the solid line and the experiment results by the dot line.

All the results described earlier show that the LDA method is suitable for investigating the geometrical structure of the hexagonal layered crystal MoS₂. Thus, in the following, we will apply the LDA method to study the electronic structure and elastic properties of the hexagonal layered crystal MoS₂.

3.2 Electronic Properties

As is known that band structure and density of states (DOSs) play an important role in analysing the physical properties of materials. The first step to explore the electronic properties of MoS₂ is to study its band structure. Our study shows that the MoS₂ is an indirect band gap semiconductor. The calculated band gap of hexagonal structure at 0 GPa is about 0.902 eV, suggesting that MoS₂ is a narrow band gap semiconductor. This value is higher than the reported gap of 0.79 eV by Li et al. [26] and 0.8 eV by Shidpour and Manteghian [31], but it more agrees with the value of 1.03 eV presented by Nayak et al. with the Perdew, Burke and Emzerhof - Generalized Gradient Approximation (PBE-GGA) methods [27] and the experimental value of 1.29 eV presented by Peelaers and Van de Walle [32]. The existing difference of band gap in our result may be addressed to the poor description of conventional density functional theory, the presence of artificial self-interaction, and the absence of the derivative discontinuity in the exchange-correlation potential of LDA methods.

The DOS is also important for understanding the bonding properties. Therefore, the total DOS (TDOS) and PDOS of MoS₂ are obtained at 0 and 14 GPa. As is clearly shown in Figure 3a, at 0 GPa, the lower valence bands (located at −15.7 ∼ −11.8 eV) are composed by S-3s states, and the upper valence bands (located at −7.404 ∼ 0 eV) are formed by S-3p and Mo-4d states, while the whole valence bands are mainly composed by the Mo-4d and S-3s states. The conduction bands are mainly from S-3p and Mo-4d states that located at 0.5 ∼ 9.6 eV. The width of S-3p band is 6.475 eV in valence band and 3.678 eV in conduction band. In addition, the width of the Mo-4p band is 6.64 eV in valence band and 3.908 eV in conduction band. Our results are consistent with those reported by Lebegue and Eriksson [33]. At 14 GPa (Fig. 3b), the lower valence bands are composed by S-3s states (located at −15.9 ∼ −12.8 eV) and Mo-4s states (located at −16.4 ∼ −12.6 eV); the upper valence bands are formed by S-3p (located at −8.2 ∼ −1.1 eV), Mo-4p (located at −8.4 ∼ −0.3 eV), and Mo-4d (located at −8.2 ∼ −0.5 eV) states, while the whole valence bands are mainly composed by the Mo-4p, Mo-4d, and S-3p states. The conduction bands are mainly from S-3p states (located at 0.6 ∼ 4.6 eV), Mo-4s states (located at 5.1 ∼ 8.9 eV), Mo-4p states (located at 0.1 ∼ 4.9 and 5.4−9.2eV), and Mo-4d states (located at 0.1 ∼ 4.9 eV). The width of S-3p band is 7.045 eV in valence band and 4.026 eV in conduction band. In addition, the width of the Mo-4p band is 8.126 eV in valence band and 9.058 eV in conduction band. From the two figures, it can be seen that the Mo-4s and Mo-4p states also contribute to the valence and conduction bands, and the whole valence bands expand towards the lower-energy area with the pressure changing from 0 to 14 GPa.

Figure 3:

Calculated TDOS and PDOS of MoS₂ at (a) 0 and (b) 14 GPa.

We make further investigations on the bonding properties of MoS₂ by the Mulliken overlap population [34], the Hirshfeld analysis, and electron density map. Mulliken’s charge Q_m(A) of A atom and the overlap population n_m(AB) for AB bond are defined as follows:

(3)Qm(A) = ∑kW(k)∑νAPμν(k)Sμν(k), (3)

(4)nm(AB) = ∑kW(k)∑μA∑νBPμν(k)Sμν(k), (4)

where P_μν(k) and S_μν(k) are the density matrix and the overlap matrix, respectively. μ and ν are the orbital on A and B, respectively. W(k) is the weight associated the calculated k points in the Brillouin zone. Usually, the magnitude and sign of Q_m(A) characterise the iconicity of A atom in the crystal, and n_m(AB) can be used to measure the covalent bonding strength approximately. In Table 1, we present the calculated Mulliken and Hirshfeld charges of MoS₂. Our Mulliken population analysis is consistent with Fu et al. from the LDA method [35]. Both results imply that the Mo atom obtains electrons, while S atoms lost electrons in the MoS₂, which is obviously unreasonable. But this can certify the conclusions that Mulliken populations seem to give an unreasonable physical picture of the charge distribution in compounds and Mulliken populations can have unphysical negative values [38]. Our Hirshfeld charge analysis implies that every Mo atom lost 0.2 electrons, while S atoms obtain 0.1 electrons. Though they seriously deviate from their nominal valence, these huge differences lie in that the Hirshfeld analysis is not a strict way in deciding the valence of an atom [39]. So, the Hirshfeld analysis is a more reasonable result. Table 1 shows the calculated bond length of Mo−S bonds 2.387 Å, which is well consistent with the other theoretical values of 2.43 [36] and 2.384 Å [37].

To display these bonding properties more visually, we plot the charge density difference of MoS₂ in Figure 4 at both 0 (a) and 14 GPa (b). In general, many isolated atoms are connected together to form a solid and liquid through such as chemical bond. The electron density difference here denotes the difference in the electron density between bonded atoms and isolated atoms. The figures reveal that there is a significant electron density in the core regions of Mo atoms but a less value in interstitial areas. The electronic charge density difference map shows some charges begin to accumulate on the red regions between the site of Mo and S, which shows the covalent characteristic of Mo–S. This result agrees with the analysis of the Milliken population. In addition, it can be seen that Mo atoms locate at the blue regions, implying that Mo atoms lost electrons. The two figures show no difference about the electron density between 0 and 14 GPa, but the electron charge becomes more and more converged between two neighbouring layers at 14 GPa. The result may be due to that the existing interlayer van del Waals (vdW) interaction increases as the pressure increases.

Figure 4:

Electronic charge density difference map calculated of MoS₂ at (a) 0 and (b) 14 GPa.

3.3 Elastic Properties

Elastic stiffness tensor is related to stress tensor and strain tensor by Hooke’s law. To calculate the elastic constants under pressures, we assume the symmetry-dependent strains to be non-volume conserving, which are appropriate for the calculations of the elastic sound velocities. The elastic constants C_ijkl with respect to the finite strain variables are defined as [40]

(5)Cijkl = (∂σij(x)∂ekl)X, (5)

where σ_ij and e_kl are the applied stress and Eulerian strain tensors, and X and x are the coordinates before and after the deformation, respectively. When the pressure P is considered,

(6a)cijkl = cijkl + P2(2δijδkl − δilδjk − δikδjl), (6a)

(6b)cijkl = (1V(x)∂2E(x)∂eij∂ekl)X, (6b)

where c_ijkl denotes the second-order derivatives with respect to the infinitesimal strain. Elastic tensor C_ijkl has generally 21 independent components. However, this number is greatly reduced when taking into account the symmetry of crystal. For the hexagonal structure MoS₂, it is reduced to five components, i.e., C₁₁, C₁₂, C₁₃, C₃₃, and C₄₄.

The calculated elastic constants are listed in Table 2, together with experimental [41] and previously calculated values [26, 42–44]. We obtain overall good agreement with experiment, especially for C₁₁, C₃₃, and C₄₄, which are all determined from neutron scattering data [41]. For C₁₂, our result shows a positive number, and [26, 42] reported the same positive sign, while [43, 44] obtained a negative sign. For C₁₃, the range of values reported by the different theoretical calculations is large, but our results are in more agreement with the calculations that inclusion of vdW interactions. The elastic constants of MoS₂ as a function of pressure are shown in Figure 5. It is found that the independent elastic constants of MoS₂ increase monotonically with pressure. C₁₁ and C₃₃ vary rapidly as pressure increases, and C₁₃, C₁₂, and C₄₄ behave moderately. The obtained elastic constants are a little different from the theoretical values by Li et al. [26]. This may be due to the small difference of the selected exchange-correlation energy that the Ceperley–Alder type parameterised by Perdew and Zunger within the LDA is used in our calcaulations, while Li et al. selected the generalised gradient approximation with PBE scheme. The criteria for the mechanical stability for hexagonal structure under pressure P are [45]

Figure 5:

The calculated elastic constants C_ij of MoS₂ as a function of pressure P. The solid lines represent our results, and the dot lines represent the calculated results by Li et al. from GGA-PBE method.

(7)C˜44 > 0, C˜11 > |C˜12|, C˜33(C˜11 + C˜12) > 2C˜132, (7)

where C˜∂∂ = C∂∂ − P(∂ = 1, 3, 4) and C˜1∂ = C1∂ + P(∂ = 2, 3, 4). The calculated elastic constants C_ij of MoS₂ satisfy the criteria, suggesting that it is mechanically stable at these applied pressures.

The bulk modulus B, shear modulus G, and Young’s modulus Y can affect physical quality of materials. Bulk modulus B is closely relative to the cohesive energy or bonding energy of atoms in crystals, so it can be used as a measure to describe the average atomic bond strength [46]. Shear modulus G is an important parameter related to the hardness of a material [47]. While Young’s modulus Y is a gauge of the stiffness of a solid, which is in a function of atomic bond strength [48]. For the hexagonal structure MoS₂, the bulk modulus B and the shear modulus G can be deduced from [49]

(8a)BV = 19[2(C11 + C12) + C33 + 4C13], (8a)

(8b)GV = 130(C11 + C12 + 2C33 − 4C13 + 12C44 + 12C66), (8b)

(8c)BR = (C11 + C12)C33 − 2C132C11 + C12 + 2C33 − 4C13, (8c)

(8d)GR = 52((C11 + C12)C33 − 2C132)C44C663BVC44C66 + ((C11 + C12)C33 − 2C132)(C44 + C66), (8d)

where C₆₆=(C₁₁ − C₁₂)/2. Young’s modulus Y and Poisson’s ratio σ can be calculated by

(9)Y = 9BG3B + G, σ = 12(1 − Y3B), (9)

where G and B are expressed as: B=(B_V +B_R)/2 and G=(G_V+G_R)/2. The calculated results are listed in Table 3. The calculated bulk modulus B at 0 GPa is in agreement with experimental data [25] and other theoretical data [26]. It can be seen from Table 3 that the elastic moduli increase gradually with the increase in pressure, that is to say, increasing pressure can improve the hardness of the materials.

Based on the assumptions that bulk modulus B represents the resistance to fracture deformation and shear modulus G represents the resistance to plastic deformation, Pugh [50] proposed the ratio B/G to distinguish ductility and brittleness of a material. The critical value that separates ductile and brittle materials is about 1.75. If B/G>1.75, the material exhibits in a ductile manner; otherwise, it behaves in a brittle manner. In this work, the calculated values of B/G of MoS₂ at ambient pressure is 1.09, indicating that MoS₂ is brittle at 0 GPa. When the applied pressure increases, the values of B/G become bigger, that is to say, pressure can reduce the brittleness of MoS₂. To further investigate the stiffness of MoS₂, the microhardness parameter (H) denoting the resistance of the physical object against compression of the contacting parts is proposed [51], H = (1 − 2σ)Y6(1 + σ). The calculated microhardness of MoS₂ as a function of pressure is also presented in Table 3. It turns out that MoS₂ is a low-stiff material. With the increase in pressure, the microhardness of MoS₂ becomes large.

As is known, the Debye temperature is an important fundamental parameter and is closely related to many physical properties of solids, such as the specific heat and melting temperature. From our elastic constants, we can obtain the elastic Debye temperature (Θ). The Debye temperature can be estimated from [52]

(10)ΘD = ℏκ[3n4π(NAρM)]1/3Vm, (10)

where ℏ is Planck’s constant, κ is Boltzmann’s constant, N_A is Avogadro’s number, n is the number of atoms per formula unit, M is the molecular mass per formula unit, ρ is the density, and V_m is obtained from [50]

(11)Vm = [13(2VS3 + Vl3)]−1/3, (11)

V_S and Vl are the shear and longitudinal sound velocities, respectively. The average shear and longitudinal sound velocities can be calculated from Navier’s equations

(12)VS = Gρ, Vl = B + (4/3)Gρ (12)

At 0 K and 0 GPa, we yield Θ = 278.35 K from the elastic constants of MoS₂. The obtained Debye temperatures of MoS₂ are also listed in Table 3. It is noted that as pressure increases, the Debye temperatures increase monotonically. For materials, usually, the higher the Debye temperature, the larger the microhardness. Therefore, the Debye temperatures also show a larger microhardness with an increasing pressure. The obtained longitudinal, shear, and average wave velocities are illustrated in Figure 6. It is shown that with the increasing pressure, the longitudinal wave velocities V_l increases fast at first, and then gently. The shear wave velocity V_S and average wave velocity V_m change slowly with the elevated pressure.

Figure 6:

Wave velocities of MoS₂ as a function of pressure. Longitudinal wave velocity V_l (the red line), average wave velocity V_m (the blue line), and shear wave velocity V_S (the black line).

To quantify the mechanical anisotropy of MoS₂, we calculate the bulk modulus B_a along a axis with the c axis fixed and the bulk modulus B_c along c axis with the a axis fixed, which are defined as follows [53]

(13a)Ba = adpda = Λ2 + α, (13a)

(13b)Bc = cdpdc = Baα, (13b)

where Λ=2(C₁₁ + C₁₂) + 4C₁₃α + C₃₃α² and α = C11 + C12 − 2C13C33 − C13. The calculated B_a, B_c and ratio B_a/B_c are presented in Figure 7, which indicates the mechanical anisotropic of MoS₂. The calculated directional bulk modulus of MoS₂ suggests that their values are larger along a axis than along c axis, which indicates that the compressibility along the a axis is smaller. The ratio B_a/B_c of MoS₂ is 7.067 at 0 GPa. It has a trend of gradual decline as the pressure increases. It is shown that the mechanical behaviour of MoS₂ at zero pressure is of strong anisotropy, but when the applied pressure increases, the anisotropy gradually decreases.

Figure 7:

Bulk modulus B_a (the blue line) and B_c (the red line) as well as B_a/B_c (the black line) versus pressure along the a and c axes.