A First-Principles Study of Nonlinear Elastic Behavior and Anisotropic Electronic Properties of Two-Dimensional HfS2

We utilize first principles calculations to investigate the mechanical properties and strain-dependent electronic band structure of the hexagonal phase of two dimensional (2D) HfS2. We apply three different deformation modes within −10% to 30% range of two uniaxial (D1, D2) and one biaxial (D3) strains along x, y, and x-y directions, respectively. The harmonic regions are identified in each deformation mode. The ultimate stress for D1, D2, and D3 deformations is obtained as 0.037, 0.038 and 0.044 (eV/Ang3), respectively. Additionally, the ultimate strain for D1, D2, and D3 deformation is obtained as 17.2, 17.51, and 21.17 (eV/Ang3), respectively. In the next step, we determine the second-, third-, and fourth-order elastic constants and the electronic properties of both unstrained and strained HfS2 monolayers are investigated. Our findings reveal that the unstrained HfS2 monolayer is a semiconductor with an indirect bandgap of 1.12 eV. We then tune the bandgap of HfS2 with strain engineering. Our findings reveal how to tune and control the electronic properties of HfS2 monolayer with strain engineering, and make it a potential candidate for a wide range of applications including photovoltaics, electronics and optoelectronics.


Introduction
The rise of two-dimensional (2D) materials began in 2004 with a focus on graphene sheets by Novoselov and Geim [1]. Graphene is a 2D layer of sp 2 -bonded carbons as the first prototype of 2D layered materials, which is viewed as an ideal material for a wide range of applications including photonics, THz electronics, nonlinear optics, sensors, and transparent electrodes [2][3][4][5]. 2D materials have been intensively researched for the next generation of ultrathin and flexible electronic and optoelectronic devices, including transistors, phototransistors, solar cells, and light-emitting diodes (LEDs) [6][7][8][9]. These materials have historically been one of the most extensively studied classes of materials due to their wealth of significant physical phenomena, which can occur when charge and heat transports are confined to a 2D surface [10][11][12].
Recently, atomically thin 2D materials, such as graphene, hexagonal boron nitride (h-BN), and the transition-metal dichalcogenides (TMDs) have received a lot of interests due to their unique electronic and optoelectronic properties. The TMDs with a general formula of MX 2 (M = transition metal, X = S, Se, Te) are particularly an interesting class of 2D materials comprising both metalic and

Computational Details
We performed DFT calculations with the Quantum ESPRESSO package [28,29] using projector-augmented wave (PAW) method. We used the Perdew-Burke-Ernzerhof exchange-correlation functional, revised for solids (PBEsol) along with the projector-augmented wave (PAW) potentials for the selfconsistent total energy calculations and geometry optimization [30]. Eighteen valence electrons of Hf atoms (4f 14 5d 2 6s 2 ) and six valence electrons of S atoms (3s 2 4p 4 ) were included in the computations. For the plane-wave expansion, the cutoff energy after convergence is set to 880 eV. The Brillouin zone sampled using a 20 × 20 × 1 Monkhorst-Pack k-point grid [31]. Atomic positions were relaxed until the energy differences are converged within 10 −6 eV and the maximum Hellmann-Feynman force on any atom is below 10 −6 eV. A vacuum of 15 Å along the c direction was included to safely avoid the interaction between the periodically repeated structures. Under various deformation tensors, the total energy of the system is calculated and led to achieve energy-strain curves. Here, the strained energy per atom is defined as below: Where E tot , and E 0 are the total energy of the strained and unstrained HfS 2 monolayers, respectively. n is also the number of atoms in the unit cell. The DFT simulation calculates the true or Cauchy stresses, σ, which for the HfS 2 monolayer should be expressed as a 2D force per length with the unit of N/m by taking the product of the Cauchy stresses (with the unit of N/m 2 ) and the super-cell thickness of 15 Å. The Cauchy stresses are related to the second Piola-Kirchhoff (PK2) stresses Σ as [32]: where, F is the deformation gradient tensor [32], J is the determinant of F, and σ is the true stress with the unit of N/m. In continuum theory of elasticity [32], and finite element method [33], the second P-K stress is employed to explore the impact of large forces on the mechanical behavior of materials [34].
Polynomial fitting of the resultant second P-K strain-stress curves on the DFT results has been conducted to calculate the second-, third-, and fourth-order elastic constants using continuum theory of elasticity.
To identify the elastic constants, the obtained strain-stress curves from DFT calculations are fitted to the constitutive equations of the continuum theory of elasticity. The second-order elastic constants (C 11 , C 12 , C 12 , and C 66 ) are the representative of the linear elastic response of the structure, while the higher-order (third-, and fourth-s order) constants are essential to the study of nonlinear elastic behavior of materials, as Wei et al. [35], Peng et al. [32], and Faghihnasiri et al. [6][7][8][9] described it completely for Boron nitride, graphene, and borophene monolayers, respectively. To identify the elastic constants, suitable deformations should be selected to facilitate the calculation of these constants directly from the stress-strain curves. For this purpose, three different types of deformations modes (strain tensors (D 1 , D 2 , and D 3 ) ) is this study, which are previously defined by Wei et al. [35], Peng et al. [32], and Faghihnasiri et al. [9].
Additionally, using the second-order elastic constants, bulk modulus (K), shear modulus (G), Young's modulus (Y) and Poisson's ratio (ν) can calculated in x and y directions. We utilize the following formulas for 2D materials to calculate these parameters [36,37]:

Atomic Structure
HfS 2 monolayer is in a hexagonal phase with an inversion center at the Hf atom sites. Each Hf atom is bounded to six S atoms. The unit cell consists of one Hf atom, and two S atoms. The geometric structure of a HfS 2 monolayer is depicted in Figure 1. The obtained optimized lattice constant of HfS 2 is 3.64 Å, which is in good agreement with the lattice constant of bulk HfS 2 , which is reported as 3.63 Å using PBE based calculations [38], and 3.61 Å utilizing vdW-TS/HI method [39]. As can be seen in Table 1, the optimized lattice constant of HfS 2 is in a geed agreement with reported lattice constant of similar 2D materials including HfSe 2 , ZrS 2 , ZrSe 2 , GaS, GaSe, and InSe in the literature. The bond length between Hf and S atoms is also calculated as 2.55Å (Figure 1b). This value is in good agreement with reported data in the literature (2.59 Å) [40]. The bond angle between Hf and S atoms (atoms No. 1, 2, and 3) is 88.80 • . Nanomaterials 2020, 10, 446 3 of 13 second P-K stress is employed to explore the impact of large forces on the mechanical behavior of materials [34]. Polynomial fitting of the resultant second P-K strain-stress curves on the DFT results has been conducted to calculate the second-, third-, and fourth-order elastic constants using continuum theory of elasticity.
To identify the elastic constants, the obtained strain-stress curves from DFT calculations are fitted to the constitutive equations of the continuum theory of elasticity. The second-order elastic constants (C11, C12, C12, and C66) are the representative of the linear elastic response of the structure, while the higher-order (third-, and fourth-s order) constants are essential to the study of nonlinear elastic behavior of materials, as Wei et al. [35], Peng et al. [32], and Faghihnasiri et al. [6][7][8][9] described it completely for Boron nitride, graphene, and borophene monolayers, respectively. To identify the elastic constants, suitable deformations should be selected to facilitate the calculation of these constants directly from the stress-strain curves. For this purpose, three different types of deformations modes (strain tensors (D1, D2, and D3) ) is this study, which are previously defined by Wei et al. [35], Peng et al. [32], and Faghihnasiri et al. [9].
Additionally, using the second-order elastic constants, bulk modulus (K), shear modulus (G), Young's modulus (Y) and Poisson's ratio (ν) can calculated in x and y directions. We utilize the following formulas for 2D materials to calculate these parameters [36,37]:

Atomic Structure
HfS2 monolayer is in a hexagonal phase with an inversion center at the Hf atom sites. Each Hf atom is bounded to six S atoms. The unit cell consists of one Hf atom, and two S atoms. The geometric structure of a HfS2 monolayer is depicted in Figure 1. The obtained optimized lattice constant of HfS2 is 3.64 Å, which is in good agreement with the lattice constant of bulk HfS2, which is reported as 3.63 Å using PBE based calculations [38], and 3.61 Å utilizing vdW-TS/HI method [39]. As can be seen in Table 1, the optimized lattice constant of HfS2 is in a geed agreement with reported lattice constant of similar 2D materials including HfSe2, ZrS2, ZrSe2, GaS, GaSe, and InSe in the literature. The bond length between Hf and S atoms is also calculated as 2.55Å (Figure 1b). This value is in good agreement with reported data in the literature (2.59 Å) [40]. The bond angle between Hf and S atoms (atoms No. 1, 2, and 3) is 88.80°.   PBEsol (This work) 3.64 GGA [40] 3.54 HSE [40] 3.53 LDA [40] 3.38 PBE (bulk) [41] 3.54 PBE (bulk) [38] 3.63 vdW-TS/HI [39] 3

Mechanical Properties
The energy-strain curves for HfS 2 monolayer under three types of deformations, namely uniaxial strain along x (D 1 ), y (D 2 ) and biaxial strain along x-y (D 3 ), is analyzed and demonstrated in Figure 2. It is evident that E s differs from the applied strain along with x and y directions. For the tensile and compressive strains through all three modes, E s becomes asymmetric. The strained energy is a quadratic function of strain in the range between −3% ≤ η ≤ 5% for the uniaxial strain along x direction. For the uniaxial strain along y direction and biaxial strain, these harmonic region ranges between −7% ≤ η ≤ 3% and −2% ≤ η ≤ 2%, respectively.

Material
Method Lattice constant

Mechanical Properties
The energy-strain curves for HfS2 monolayer under three types of deformations, namely uniaxial strain along x (D1), y (D2) and biaxial strain along x-y (D3), is analyzed and demonstrated in Figure 2. It is evident that differs from the applied strain along with x and y directions. For the tensile and compressive strains through all three modes, becomes asymmetric. The strained energy is a quadratic function of strain in the range between −3% ≤ ≤ 5% for the uniaxial strain along x direction. For the uniaxial strain along y direction and biaxial strain, these harmonic region ranges between −7% ≤ ≤ 3% and −2% ≤ ≤ 2%, respectively.
In all three deformation modes, the total energy of the system increases with increasing the applying strain. As can be seen in Figure 2, The changes of Es with strain for D1 and D2 deformation modes is almost the same. However, for D3 deformation mode, the rate of energy change with strain is much higher (Figure 2).  Figure 3 shows the strain-stress relations obtained from the DFT calculations as well as the fits to these by the equations of continuum theory of elasticity. As can be seen in this figure, the stressstrain curves are depicted for all D1, D2, and D3 deformation modes. In all three deformation modes, the total energy of the system increases with increasing the applying strain. As can be seen in Figure 2, The changes of E s with strain for D 1 and D 2 deformation modes is almost the same. However, for D 3 deformation mode, the rate of energy change with strain is much higher (Figure 2). Figure 3 shows the strain-stress relations obtained from the DFT calculations as well as the fits to these by the equations of continuum theory of elasticity. As can be seen in this figure, the stress-strain curves are depicted for all D 1 , D 2 , and D 3 deformation modes. The maximum value in the stress-strain curves shows the ultimate tensile strength (Σm) of the material, in which a material can suffer the maximum respective ηm without damaging ( Figure 4). The ultimate strain reflects the intrinsic bonding strengths and acts as a lower limit of the critical strain. Additionally, the values obtained for the ultimate stress and ultimate strain is given in Table  2. Beyond the ultimate strain, the materials will get in a metastable state, which ends up with fracture [42]. The DFT results for strains below the ultimate strain are used to determine the higher-order The ultimate strain reflects the intrinsic bonding strengths and acts as a lower limit of the critical strain. Additionally, the values obtained for the ultimate stress and ultimate strain is given in Table 2. Beyond the ultimate strain, the materials will get in a metastable state, which ends up with fracture [42]. The DFT results for strains below the ultimate strain are used to determine the higher-order elastic constants. As can be seen in Figure 4, the stress increases linearly with increasing strain, within the harmonic (elastic) region. Under larger strains, for the prediction of strain-stress curves, the system is in the anharmonic region in which higher-order terms must be perceived as well. As mentioned previously, the system transits from elastic to the plastic region, when exposed to higher strains. Eventually, Table 3 prepares the nonzero second-, third-and fourth-order elastic constants (SOEC, TOEC and FOEC, respectively) for the HfS2 monolayer.  To comprehensively understand the magnitudes of elastic constants obtained in this work for HfS2, Table 4 presents multiple comparisons between our findings and the other reported elastic constants for similar structures. Furthermore, 2D Young's moduli (in-plane stiffness) along x and y directions ( , ), Poisson's ratio along x and y directions ( , ), 2D shear modulus (G 2D ), and the 2D bulk modulus (K), are tabulated in Table 5 and validated by those reported for other structurally similar compounds in the literature.  As can be seen in Figure 4, the stress increases linearly with increasing strain, within the harmonic (elastic) region. Under larger strains, for the prediction of strain-stress curves, the system is in the anharmonic region in which higher-order terms must be perceived as well. As mentioned previously, the system transits from elastic to the plastic region, when exposed to higher strains. Eventually, Table 3 prepares the nonzero second-, third-and fourth-order elastic constants (SOEC, TOEC and FOEC, respectively) for the HfS 2 monolayer. To comprehensively understand the magnitudes of elastic constants obtained in this work for HfS 2 , Table 4 presents multiple comparisons between our findings and the other reported elastic constants Nanomaterials 2020, 10, 446 7 of 13 for similar structures. Furthermore, 2D Young's moduli (in-plane stiffness) along x and y directions (Y 2D

Strain-Stress Relationship
x , Y 2D y ), Poisson's ratio along x and y directions (v 2D x , v 2D y ), 2D shear modulus (G 2D ), and the 2D bulk modulus (K), are tabulated in Table 5 and validated by those reported for other structurally similar compounds in the literature. Table 4. Elastic constants C 11 and C 12 obtained in this work along with those reported in a previous works (in units of N/m).  [39] 131. 47 25.63 ZrSe 2 [39] 104.62 21.31 HfS 2 [39] 141.98 25.95 HfSe 2 [39] 116.88 22.30 Ref [39] is in the bulk structure.

Electronic Properties
First, we studied the electronic properties of HfS 2 in the absence of strain. Figure 5 shows the band structure of the strained structure of HfS 2 , which is obtained from PBEsol calculations. As can be seen in Figure 5, HfS 2 is a semiconductor with an indirect bandgap of 1.12 eV. This value is compared with other methods of HfS 2 (E LDA g = 1.07 eV [40], E GGA g = 1.15 eV [43], E HSE06 g = 2.02 eV [44], and E GW g = 2.45 eV [45]), where it is in a good agreement with E GGA g = 1.15 eV [43]. As can be seen in Table 4, The predicted gap energies by GW and HSE06 are larger than the predicted gap energy by PBE approximation. Since the HSE06 functional can accurately predict enthalpies of formation, ionization potentials, and electron affinities for lattice constants and band gaps of solids in general, the predicted gap energy by HSE06 is larger and more accurate than the predicted gap energy by PBE [46]. Also, in GW calculations, the excitonic effects are included and due to the fact that the excitonic effects are significant in 2D semiconductors, the predicted gap energy by GW is twice larger than the predicted gap energy by GGA approximation for the HfS 2 monolayer (Table 4) [47].
In our predicted electronic structure ( Figure 5), it can be seen that the conduction-band minimum (CBM) and valence-band maximum (VBM) are located at M and Γ points, respectively. Additionally, we examine the projected density of states (PDOS), which demonstrates the contribution of orbitals in the valence and conduction bands of the material. In the conduction band, Hf-5d orbitals have the maximum contribution. In contrast, the 3p orbital of S atom makes a greater contribution in the valance band near the Fermi level, as shown in Figure 5b.
In the next step, we investigated the electronic behavior of HfS 2 under different strains. The bandgap variation of HfS 2 monolayer under these three types of strain are shown in Figure 6a. as can be seen in this figure, for all deformation modes, the bandgap decreases when the compressive strain increases 0% to 10%, while it increases when the tensile strain increases from 0% to 10%. In this strain range, we note that the rate of the gap energy decrease/increase with strain is almost the same for D1 and D2 deformation modes. From 10% to 16% tensile strain, the energy gap decreases for all deformation modes. For the strain range from 16% to 20%, the energy gap increases by increasing strain for D1 and D2 deformations, while it keeps decreasing for D3. For the tensile strain ranging from 20% to 30%, the gap energy is almost constant with strain for D1 and D2, while it keeps decreasing for D3. At 22%, 24%, 28% strains along x direction (D 1 ), and at 12%, 14%, 16%, 20%, 22% strains along y direction (D 2 ), the bandgap becomes direct. Under biaxial strain (D 3 ), the semiconducting monolayer maintains its indirect nature, while its energy gap increases with increasing strain and reaches it maximum value (1.79 eV) at 10% strain. For the uniaxial strain in both D 1 and D 2 cases, when the compressive strain increases from zero to 2%, the energy gap decreases and then gradually increases by increasing the amount of the compressive strain from 2% to 5%. Then, by increasing the strain from 5% to 10%, the energy gap reduces again. For D3 case, the energy gap decreases by increasing the compressive strain from 0 to 10% deformation. As can be seen in Figure 6a, the material becomes metal (E g = 0 eV) at 10% compressive strain in all three deformation modes. The value of energy gap is the same (1.6 eV) for D 1 and D 3 at 18% strain, and for D 2 at 17.5% strain. As can be observed in Figure 4, at 18% D1 and 17.5% D2 types of deformation, the material is in plastic region, while at 18% D 3 type deformation, the material shows elastic deformation (Figure 6b,c,d). As can be seen in Figure 6d, the elastic to plastic transition for D3 occurs at 22% strain. In the next step, we investigated the electronic behavior of HfS2 under different strains. The bandgap variation of HfS2 monolayer under these three types of strain are shown in Figure 6a. as can be seen in this figure, for all deformation modes, the bandgap decreases when the compressive strain increases 0% to 10%, while it increases when the tensile strain increases from 0% to 10%. In this strain range, we note that the rate of the gap energy decrease/increase with strain is almost the same for D1 and D2 deformation modes. From 10% to 16% tensile strain, the energy gap decreases for all deformation modes. For the strain range from 16% to 20%, the energy gap increases by increasing strain for D1 and D2 deformations, while it keeps decreasing for D3. For the tensile strain ranging from 20% to 30%, the gap energy is almost constant with strain for D1 and D2, while it keeps decreasing for D3. At 22%, 24%, 28% strains along x direction (D1), and at 12%, 14%, 16%, 20%, 22% strains along y direction (D2), the bandgap becomes direct. Under biaxial strain (D3), the semiconducting monolayer maintains its indirect nature, while its energy gap increases with increasing strain and reaches it maximum value (1.79 eV) at 10% strain. For the uniaxial strain in both D1 and D2 cases, when the compressive strain increases from zero to 2%, the energy gap decreases and then gradually increases by increasing the amount of the compressive strain from 2% to 5%. Then, by increasing the strain from 5% to 10%, the energy gap reduces again. For D3 case, the energy gap decreases by increasing the compressive strain from 0 to 10% deformation. As can be seen in Figure 6a, the material becomes metal (Eg = 0 eV) at 10% compressive strain in all three deformation modes. The value of energy gap is the same (1.6 eV) for D1 and D3 at 18% strain, and for D2 at 17.5% strain. As can be observed in Figure 4, at 18% D1 and 17.5% D2 types of deformation, the material is In the next step, we investigated the electronic behavior of HfS2 under different strains. The bandgap variation of HfS2 monolayer under these three types of strain are shown in Figure 6a. as can be seen in this figure, for all deformation modes, the bandgap decreases when the compressive strain increases 0% to 10%, while it increases when the tensile strain increases from 0% to 10%. In this strain range, we note that the rate of the gap energy decrease/increase with strain is almost the same for D1 and D2 deformation modes. From 10% to 16% tensile strain, the energy gap decreases for all deformation modes. For the strain range from 16% to 20%, the energy gap increases by increasing strain for D1 and D2 deformations, while it keeps decreasing for D3. For the tensile strain ranging from 20% to 30%, the gap energy is almost constant with strain for D1 and D2, while it keeps decreasing for D3. At 22%, 24%, 28% strains along x direction (D1), and at 12%, 14%, 16%, 20%, 22% strains along y direction (D2), the bandgap becomes direct. Under biaxial strain (D3), the semiconducting monolayer maintains its indirect nature, while its energy gap increases with increasing strain and reaches it maximum value (1.79 eV) at 10% strain. For the uniaxial strain in both D1 and D2 cases, when the compressive strain increases from zero to 2%, the energy gap decreases and then gradually increases by increasing the amount of the compressive strain from 2% to 5%. Then, by increasing the strain from 5% to 10%, the energy gap reduces again. For D3 case, the energy gap decreases by increasing the compressive strain from 0 to 10% deformation. As can be seen in Figure 6a, the material becomes metal (Eg = 0 eV) at 10% compressive strain in all three deformation modes. The value of energy gap is the same (1.6 eV) for D1 and D3 at 18% strain, and for D2 at 17.5% strain. As can be observed in Figure 4, at 18% D1 and 17.5% D2 types of deformation, the material is in plastic region, while at 18% D3 type deformation, the material shows elastic deformation ( Figure  6b,c,d). As can be seen in Figure 6d, the elastic to plastic transition for D3 occurs at 22% strain. Our findings indicate that the electronic properties of HfS2 can be effectively tuned by applying planar forces to HfS2 in different directions. Figure 7 shows the band structures of the HfS2 monolayer under uniaxial compressive and tensile strains along the x direction (D1) within the range of −10% to 30%. As can be seen in Figure 7, when the strain ranges from −10% to −8%, bands of energies cross the Fermi level so that HfS2 at these strains under D1 deformation shows a metallic behavior. By reducing the value of strain to −6%, the HfS2 monolayer becomes an indirect semiconductor with a bandgap of 0.35 eV. At the strains above 22%, the bandgap becomes direct and energy of bandgap increases up to 1.8 eV for the D1 case.
In Figure 7, the conduction-band minimum (CBM) is located at S point, while the valence-band maximum (VBM) is located at the Γ point, which is shifted to the S point by increasing the strain. These points indicate the transformation of bandgap from indirect to direct with strain engineering. For the deformed HfS2 under the tensile strain along y direction (D2), by increasing the strains from 0% to 30%, the bandgap increases from 1.12 eV to its maximum value of 2.11 eV continuously. The band structure of strain HfS2 monolayer under uniaxial compressive and tensile strain along y direction (D2) is shown in Figure 8. Under the compressive strain of D2 deformation, there is no bandgap near the Fermi level so the system is metallic. In addition, it remains as an indirect Our findings indicate that the electronic properties of HfS 2 can be effectively tuned by applying planar forces to HfS 2 in different directions. Figure 7 shows the band structures of the HfS 2 monolayer under uniaxial compressive and tensile strains along the x direction (D 1 ) within the range of −10% to 30%. As can be seen in Figure 7, when the strain ranges from −10% to −8%, bands of energies cross the Fermi level so that HfS 2 at these strains under D 1 deformation shows a metallic behavior. By reducing the value of strain to −6%, the HfS 2 monolayer becomes an indirect semiconductor with a bandgap of 0.35 eV. At the strains above 22%, the bandgap becomes direct and energy of bandgap increases up to 1.8 eV for the D 1 case. Our findings indicate that the electronic properties of HfS2 can be effectively tuned by applying planar forces to HfS2 in different directions. Figure 7 shows the band structures of the HfS2 monolayer under uniaxial compressive and tensile strains along the x direction (D1) within the range of −10% to 30%. As can be seen in Figure 7, when the strain ranges from −10% to −8%, bands of energies cross the Fermi level so that HfS2 at these strains under D1 deformation shows a metallic behavior. By reducing the value of strain to −6%, the HfS2 monolayer becomes an indirect semiconductor with a bandgap of 0.35 eV. At the strains above 22%, the bandgap becomes direct and energy of bandgap increases up to 1.8 eV for the D1 case.
In Figure 7, the conduction-band minimum (CBM) is located at S point, while the valence-band maximum (VBM) is located at the Γ point, which is shifted to the S point by increasing the strain. These points indicate the transformation of bandgap from indirect to direct with strain engineering. For the deformed HfS2 under the tensile strain along y direction (D2), by increasing the strains from 0% to 30%, the bandgap increases from 1.12 eV to its maximum value of 2.11 eV continuously. The band structure of strain HfS2 monolayer under uniaxial compressive and tensile strain along y direction (D2) is shown in Figure 8. Under the compressive strain of D2 deformation, there is no bandgap near the Fermi level so the system is metallic. In addition, it remains as an indirect In Figure 7, the conduction-band minimum (CBM) is located at S point, while the valence-band maximum (VBM) is located at the Γ point, which is shifted to the S point by increasing the strain. These points indicate the transformation of bandgap from indirect to direct with strain engineering.
For the deformed HfS 2 under the tensile strain along y direction (D 2 ), by increasing the strains from 0% to 30%, the bandgap increases from 1.12 eV to its maximum value of 2.11 eV continuously. The band structure of strain HfS 2 monolayer under uniaxial compressive and tensile strain along y direction (D 2 ) is shown in Figure 8. Under the compressive strain of D 2 deformation, there is no bandgap near the Fermi level so the system is metallic. In addition, it remains as an indirect semiconductor over the entire applied compressive strain domain when it is strained along y direction (D 2 ). However, under tensile straining, at 12% to 22% strains, the gap is direct and energy bandgap changes from 1.33 to 2.06 eV. When the compressive strain is applied (D2), the located CBM at S point moves to Γ point, and the VBM at the Γ point moves to M point. By applying tensile strain, CBM gets away from the Fermi level and the bandgap increases from 1.12 eV (at unstrained condition) to 2.11 eV (strained with 30% tensile strain). Nanomaterials 2020, 10, 446 10 of 13 semiconductor over the entire applied compressive strain domain when it is strained along y direction (D2). However, under tensile straining, at 12% to 22% strains, the gap is direct and energy bandgap changes from 1.33 to 2.06 eV. When the compressive strain is applied (D2), the located CBM at S point moves to Γ point, and the VBM at the Γ point moves to M point. By applying tensile strain, CBM gets away from the Fermi level and the bandgap increases from 1.12 eV (at unstrained condition) to 2.11 eV (strained with 30% tensile strain). In Figure 9 the band structure of strained HfS2 monolayer under biaxial strain along x-y (D3) ranging from −10% to 30% strain is shown. As can be seen in Figure 9, under compressive strain ranging from −10% to −6% strains, the bands cross the Fermi level and the system shows metallic behavior. At −6% strain and beyond (−6 < strain < 0), the HfS2 monolayer becomes an indirect semiconductor with a bandgap of 0.16 eV in −6% strain. In the tensile deformation domain, the gap energy increases increasing the stain. At 0%, 10%, 20%, and 30% strains, the bandgap becomes 1.12 eV, 1.79 eV, 1.59 eV and 1.45 eV, respectively. The CBM is located at M, Г, Г and K point for 0%, 10%, 20%, and 30% strain, respectively. Also, the VBM at 0%, 10%, 20% strains are located between M and Г points and at 30% strain, it is located between Г and K points. Graphene has a plate structure and its symmetry is not broken during deformation. Therefore, not only graphene does not suffer any bulking during deformation, but also only monotonic changes can occur on its electron properties during the biaxial straining ( Figure 10) [48]. While, the HfS2 In Figure 9 the band structure of strained HfS 2 monolayer under biaxial strain along x-y (D 3 ) ranging from −10% to 30% strain is shown. As can be seen in Figure 9, under compressive strain ranging from −10% to −6% strains, the bands cross the Fermi level and the system shows metallic behavior. At −6% strain and beyond (−6 < strain < 0), the HfS 2 monolayer becomes an indirect semiconductor with a bandgap of 0.16 eV in −6% strain. In the tensile deformation domain, the gap energy increases increasing the stain. At 0%, 10%, 20%, and 30% strains, the bandgap becomes 1.12 eV, 1.79 eV, 1.59 eV and 1.45 eV, respectively. The CBM is located at M, Г, Г and K point for 0%, 10%, 20%, and 30% strain, respectively. Also, the VBM at 0%, 10%, 20% strains are located between M and Г points and at 30% strain, it is located between Г and K points. semiconductor over the entire applied compressive strain domain when it is strained along y direction (D2). However, under tensile straining, at 12% to 22% strains, the gap is direct and energy bandgap changes from 1.33 to 2.06 eV. When the compressive strain is applied (D2), the located CBM at S point moves to Γ point, and the VBM at the Γ point moves to M point. By applying tensile strain, CBM gets away from the Fermi level and the bandgap increases from 1.12 eV (at unstrained condition) to 2.11 eV (strained with 30% tensile strain). In Figure 9 the band structure of strained HfS2 monolayer under biaxial strain along x-y (D3) ranging from −10% to 30% strain is shown. As can be seen in Figure 9, under compressive strain ranging from −10% to −6% strains, the bands cross the Fermi level and the system shows metallic behavior. At −6% strain and beyond (−6 < strain < 0), the HfS2 monolayer becomes an indirect semiconductor with a bandgap of 0.16 eV in −6% strain. In the tensile deformation domain, the gap energy increases increasing the stain. At 0%, 10%, 20%, and 30% strains, the bandgap becomes 1.12 eV, 1.79 eV, 1.59 eV and 1.45 eV, respectively. The CBM is located at M, Г, Г and K point for 0%, 10%, 20%, and 30% strain, respectively. Also, the VBM at 0%, 10%, 20% strains are located between M and Г points and at 30% strain, it is located between Г and K points. Graphene has a plate structure and its symmetry is not broken during deformation. Therefore, not only graphene does not suffer any bulking during deformation, but also only monotonic changes can occur on its electron properties during the biaxial straining ( Figure 10) [48]. While, the HfS2 Graphene has a plate structure and its symmetry is not broken during deformation. Therefore, not only graphene does not suffer any bulking during deformation, but also only monotonic changes can occur on its electron properties during the biaxial straining ( Figure 10) [48]. While, the HfS 2 structure is not a plate like structure, its electronic properties are affected by the occurred buckling during the deformation. It is expected to observe many changes in the symmetry of the structure and the distance and angle of the atoms in the strained HfS 2 monolayers. Such structural changes can greatly affect the electronic properties and therefore, the bandgap changes of HFS 2 will not be monotonic under biaxial strains ( Figure 6). Nanomaterials 2020, 10, 446 11 of 13 structure is not a plate like structure, its electronic properties are affected by the occurred buckling during the deformation. It is expected to observe many changes in the symmetry of the structure and the distance and angle of the atoms in the strained HfS2 monolayers. Such structural changes can greatly affect the electronic properties and therefore, the bandgap changes of HFS2 will not be monotonic under biaxial strains ( Figure 6).

Conclusions
In summary, the mechanical and electronic properties of the HfS2 monolayer under two uniaxial (D1 and D2) and one biaxial (D3) DFT calculations is investigated. We determined the harmonic regions in the different deformation modes. This harmonic region ranges in −3% ≤ ≤ 5% , −7% ≤ ≤ 3% and −2% ≤ ≤ 2% for D1, D2, and D3, respectively. Our findings reveal that the ultimate stress of the HfS2 monolayer for D1, D2, and D3 is 0.037 Å , 0.038 Å ,and 0.044 Å , respectively. The obtained ultimate strain is 17%, 17.5% and 21% strain for D1, D2, and D3 respectively. The high order of elastic constants including second-, third-, and fourth-order constants are calculated. The values of 2D Young's moduli along x and y directions are predicted as 83.01 N/m and 83.57 N/m, respectively. The value of Poisson's ratio along x and y directions is the same (0.17) for both D1 and D2. Moreover, the electronic properties of HfS2 show that it is a semiconductor with an indirect bandgap of 1.12 eV. The projected density of states (PDOS) indicates the conduction band, Hf-5d orbital possesses the maximum contribution, while the 3p orbital of S atom have greatest contribution in the valance band. The variation of band structure and bandgap of the HfS2 monolayer under D1, D2, and D3 deformation modes in the range of −10% to 30% are also investigated. We tuned the bandgap state (direct vs. indirect), gap energy (opening vs. shrinking), and phase transition (semiconductor-metal) by strain engineering under different deformation modes. Our findings reveal how to utilize strain engineering to make HfS2 monolayer as a suitable candidate for a wide range of applications including flexible solar cells, electronics and optoelectronics.

Conclusions
In summary, the mechanical and electronic properties of the HfS 2 monolayer under two uniaxial (D 1 and D 2 ) and one biaxial (D 3 ) DFT calculations is investigated. We determined the harmonic regions in the different deformation modes. This harmonic region ranges in −3% ≤ η ≤ 5%, −7% ≤ η ≤ 3% and −2% ≤ η ≤ 2% for D 1 , D 2 , and D 3 , respectively. Our findings reveal that the ultimate stress of the HfS 2 monolayer for D 1 , D 2 , and D 3 is 0.037 eV A 3 , 0.038 eV A 3 ,and 0.044 eV A 3 , respectively. The obtained ultimate strain is 17%, 17.5% and 21% strain for D 1 , D 2 , and D 3 respectively. The high order of elastic constants including second-, third-, and fourth-order constants are calculated. The values of 2D Young's moduli along x and y directions are predicted as 83.01 N/m and 83.57 N/m, respectively. The value of Poisson's ratio along x and y directions is the same (0.17) for both D1 and D2.
Moreover, the electronic properties of HfS 2 show that it is a semiconductor with an indirect bandgap of 1.12 eV. The projected density of states (PDOS) indicates the conduction band, Hf-5d orbital possesses the maximum contribution, while the 3p orbital of S atom have greatest contribution in the valance band. The variation of band structure and bandgap of the HfS 2 monolayer under D 1 , D 2 , and D 3 deformation modes in the range of −10% to 30% are also investigated. We tuned the bandgap state (direct vs. indirect), gap energy (opening vs. shrinking), and phase transition (semiconductormetal) by strain engineering under different deformation modes. Our findings reveal how to utilize strain engineering to make HfS 2 monolayer as a suitable candidate for a wide range of applications including flexible solar cells, electronics and optoelectronics.