Strain-induced effects on the electronic properties of 2D materials

Thanks to the ultrahigh flexibility of 2D materials and to their extreme sensitivity to applied strain, there is currently a strong interest in studying and understanding how their electronic properties can be modulated by applying a uniform or nonuniform strain. In this work, using density functional theory (DFT) calculations, we discuss how uniform biaxial strain affects the electronic properties, such as ionization potential, electron affinity, electronic gap, and work function, of different classes of 2D materials from X-enes to nitrides and transition metal dichalcogenides. The analysis of the states in terms of atomic orbitals allows to explain the observed trends and to highlight similarities and differences among the various materials. Moreover, the role of many-body effects on the predicted electronic properties is discussed in one of the studied systems. We show that the trends with strain, calculated at the GW level of approximation, are qualitatively similar to the DFT ones solely when there is no change in the character of the valence and conduction states near the gap.


Introduction
Much of the emphases on 2D materials 1-3 was born with the discovery of graphene, for which the Nobel Prize in physics was awarded to Novoselov and Geim in 2010. 4 Graphene is a 2D crystal made up of carbon atoms arranged in a hexagonal honeycomb form; it is one million times thinner than paper, almost transparent, and, at the same time, is the strongest material in the world. 5,6 Its electronic structure can be easily derived from a simple tight-binding model, which explains the presence of bands with conical dispersion intersecting at the Fermi level, thus making graphene a semimetal. 7 Massless Dirac fermions move in graphene as fast as v F ¼ 10 6 m=s at the Fermi level, and twisted bilayers of graphene have even been shown to be superconductors for very small twisting angles. 8 This peculiar linear dispersion of low-energy carriers, which can be mapped to an effective 2D Dirac Hamiltonian, is very different from the usual parabolic dispersion of bulk semiconductors. This has stimulated a lot of work to theoretically predict and experimentally observe novel physical effects in this 2D material. Nevertheless, despite the enormous interest both at the fundamental and applicative level, the lack of an electronic gap limits graphene use in applications like digital electronics, field-effect transistors, and optoelectronics at visible frequencies. For this reason, the last years have witnessed growing efforts to find a way to open its gap, leaving unmodified its peculiar electronic and mechanical properties, but also to grow and characterize other novel metallic and semiconducting two-dimensional materials beyond graphene.
Following the route of graphene, broad families of 2D materials are hence continuously developed and studied in view of their interesting physical properties and of a large number of their envisaged device-oriented applications. 9,10 Because of their atomic-scale thickness, they are characterized by weak dielectric screening, strong light-matter interaction, and highly bound excitons. Moreover, atomic and molecular doping, [11][12][13] external fields, 14 and also strain [15][16][17] may have a very deep impact on their electronic and optical properties. In particular, strain engineering is very exciting since, differently from 3-D traditional materials, 2D materials can endure remarkably large mechanical strain (up to 10%), hence creating opportunities to modulate their physical properties for interesting device applications.
Silicene, [18][19][20] the silicon-based counterpart of graphene, represents the first exciting material merging the exceptional physical properties of graphene with the simplicity of easily integrating it in already existing and largely developed silicon-based technology. 21,20 A field-effect transistor has been reported at room temperature. 22 Further interest in silicene arises from its predicted nontrivial topological properties. 23 Freestanding ideal silicene presents a buckled honeycomb structure. 19 Like graphene, silicene has a semimetallic behavior and possesses (in the absence of spinorbit corrections) massless fermions at the k-point of the brillouin zone (BZ). 24,25 Other members of the so-called X-enes family (borophene, germanene, stanene, phosphorene, arsenene, antimonene, bismuthene, and tellurene) are also of particular interest for their excellent physical, chemical, electronic, and optical properties. 26 The evidence that bulk group-III nitrides are among the most important materials for solid-state lighting, as witnessed by the Nobel prize awarded in 2014 to Akasaki, Amano, and Nakamura, stimulated in the last years several theoretical 27,28 and experimental 29-31 studies on 2D honeycomb III-N sheets. From the experimental side, their growth is very challenging because, similar to the case of silicene, no simple route to mechanical or chemical exfoliation can be used, due to their 3-D wurtzite structures (only BN crystallizes in the hexagonal layered form in bulk). However, several promising experimental attempts to realize 2D III-nitrides have been reported so far.
Among layered materials, for which exfoliation is possible, transition metal dichalcogenide (TMD) occupies a prominent place in recent worldwide research. [32][33][34][35] They are of type MX 2 , where one layer of transition metal atoms (M) is sandwiched between two layers of chalcogen (X) atoms, crystallizing mainly in the hexagonal or rhombohedral forms with metal atoms having octahedral or trigonal prismatic coordination. Among this broad family, group-VI TMD (MX 2 with M ¼ Mo, W, and X ¼ S, selenium, tellurium (Te)) has received a lot of attention due to their tunable bandgap, strong light-matter interaction, 36,37 large spinorbit effects, strong influence of doping, functionalization, external field, and strain. In the hexagonal form, they are direct gap semiconductors when monolayer (ML) while exhibiting indirect gaps when BL and thicker multilayers, being promising materials for flexible electronics, light emission, energy storage, solar energy conversion, as well as electrochemical catalysis and biosensors. [38][39][40][41] In this article, we use first-principle DFT calculations to study the role of biaxial uniform strain on the structural and electronic properties of 2D materials beyond graphene: silicene, group III-nitrides, and TMD. In particular, we focus on the behavior of their band structure, electronic gap, workfunction (WF), ionization potential (IP), and electron affinity (EA) and perform a comparative analysis for nonpolar (graphene and silicene), polar (aluminum nitride (AlN) and gallium nitride (GaN)) honeycomb lattices, and TMD (molybdenum disulfide (MoS 2 ) and molybdenum ditelluride (MoTe 2 )) in the monolayer and bilayer form.

Methods and computational details
All the DFT calculations to obtain the structural and electronic properties have been performed using the Quantum-Espresso code 42 within the local density approximation (LDA), 43 perdew becke ernzerhrof (PBE), 44 XC functional for nitrides (TMD). For the case of the two Xenes (graphene and silicene), we performed calculations using both LDA and PBE with the specific goal to look at the effect of local or semilocal XC on the calculated values of the work functions. The van der waals (vdW) correction 45 has been applied to take into account the interaction between the layers in TMD BL.
The self-consistent density for graphene and silicene was computed using a 15 Â 15 Â 1 k-point sampling for the ground state with plane wave cutoff of 100 Ry. For nitrides, we used a nonshifted 18 Â 18 Â 1 k-point sampling for the ground state and a plane wave cutoff of 100 (200) Ry for the structural optimization and band structure of AlN (GaN). For TMD, we sampled the BZ using a uniform 18 Â 18 Â 1 mesh and used a plane wave cutoff of 120 Ry.
To simulate isolated layers, vacuum sizes of 15 and 20 Å have been chosen for X-enes and nitrides, respectively, while 17.8 (16.5) Å has been used for TMD-ML (BL). The geometry of each system has been relaxed at each value of strain until the forces on the atoms were less than 10 À2 eV A À1 . The values of applied strain are reported in the form e ¼ aÀa 0 a 0 , where a 0 is the equilibrium lattice constant. Positive values correspond to tensile strain, while negative values correspond to compressive strain. The equilibrium geometries at zero strain are in very good agreement with the existing literature 47,48 and the main equilibrium structural parameters are reported in Table 1.
To assess the role of many-body effects, for one of the studied systems, that is, the MoTe 2 monolayer, we went beyond the single-particle approach. On top of the DFT simulations, we performed one-shot perturbative GW calculations at 2%; 0; À2% strain using the Yambo code. 49,50 A cutoff of 240 Ry has been used for the exchange part of the quasiparticle corrections < S x À V xc >, while 8 Ry and 100 bands were employed for the correlation part of the selfenergy < S c >. To speed up the convergence with respect to the empty states, we adopted the technique described by Bruneval F and Gonze. 51 A plasmon-pole approximation for the inverse dielectric matrix has been applied. 52 Moreover, to guarantee the simulation of isolated layers, a cutoff in the bare Coulomb potential has also been used. 50,53 The k-point sampling was selected to be 42 Â 42 Â 1 in the BZ.

Results
The main goal of this work is to describe how the electronic properties of different classes of metallic and semiconducting 2D materials can be tuned via the application of a uniform biaxial strain with the aim to highlight similarities and differences among the different systems.
For all the 2D materials, we have calculated the WF obtained as the energy difference between the Fermi energy E F and the vacuum potential E vac . For 2D semiconductors, we report the electronic band structures and bandgap values at several biaxial uniform compressive and tensile strains, the IP ( IP ¼ E vac À E VBM ) and the EA ( EA ¼ E vac À E CBM ) to show their specific dependence on strain and to understand which are the different effects leading to the peculiar behaviors observed. For nitrides (TMD ML and their homobilayers), the DFT band structures have been calculated for applied strains ranging from À10% to 10% (À7% to 7% for MoS 2 ,À10% to 10% for MoTe 2 ) but, here, we report only those at zero strain and À4%, 4% (À5%, 5% for MoS 2 and MoTe 2 ). For each atomic structure at a given strain, all the quantities are determined self-consistently in a ground-state electronic structure calculation. The vacuum level is obtained as the asymptotic value of the electrostatic potential in the direction perpendicular to the layer at a far distance, 10 Å or more, from the system.

Graphene and silicene
For these two X-enes, we observe that for increasing compressive strain, the WF decreases, whereas it increases for increasing tensile strain (see Figure 1). For graphene, this behavior is in substantial agreement with the existing literature 54,55 and it has been explained in terms of a straininduced enhancement of the density of states (DOS) close to the Fermi level. 56,57 Results for silicene are qualitatively similar to the ones obtained for graphene. In the tensile region, the WF variation due to a 10% of strain is found to be 0:15 eV, in good agreement with the results of Qin et al., 58,59 where a change of about 0:2 eV has been obtained.
We observe that, for both the systems, a simple argument of plausibility of this behavior is based on the fact that as long as the materials are stretched, the interaction among the ions of the lattice decreases, thus approaching the behavior of isolated atoms in the limit of infinite tensile deformation. The ionization potentials for C and Si atoms are 11:2 and 8:1 eV, respectively, and these values are much higher than the WF of their correspondent 2D forms at equilibrium, 4:24 eV for graphene, and 4:35 eV for silicene. For this reason, it is expected that the WF of these materials, characterized by fully covalent bonds, should grow with increasing uniform tensile deformation.
To test the possible dependence of our results on the choice of the XC potential used to describe the various systems, we have analyzed the strain dependence of the WF of graphene and silicene both within the LDA and PBE scheme for the XC term. The results of this analysis are shown in Figure 1 in different colors. It is noted that the Table 1. Equilibrium lattice constant a 0 , buckling for X-enes and nitrides, and X-M-X angle for TMD.  choice of the functional affects the value of the WF by adding a constant shift of 0:24 eV for graphene and of 0:41 eV for silicene. In particular, the zero strain WF for graphene ranges from 4:24 eV in PBE to 4:48 eV in LDA, thus producing a spread of results, which is fully consistent with values in the range from 4:28 eV to 4:5 eV found in the literature. 54,55,[60][61][62] But, interestingly, we observe that the functional induced variation turns out to be a constant shift, which is independent of the strain value. Thus, we can conclude that the strain-induced variations of the WF do not depend on the choice of functional and can be assessed unambiguously within the DFT framework.

Nitrides monolayers: AlN and GaN
The 2D nitrides have the same flat honeycomb structure as graphene but with polar bonds due to the different electronegativity of the lattice atoms. To understand the trends of the IP, WF, and EA under strain, it is important to analyze their band structures. These are shown in Figure 2 for GaN and AlN at 0, À4%, and 4% uniform strain. In agreement with existing literature, 27 we have found that at zero strain, both AlN and GaN have an indirect bandgap G c À K v and the gap (increases) decreases for (compressive) tensile uniform strain, as shown in Figure 2. Interestingly, an important difference occurs between AlN and GaN. While AlN remains an indirect gap material for all the values of considered strain, GaN undergoes an indirect to direct gap transition at the compressive strain of about À4%. It is worth to mention that unstrained GaN has an indirect gap only at the DFT level, which switches to direct when the GW corrections are introduced. 27, 63 In the case of strain, the VBM switches from K to G, because the state at G goes up in energy under compressive strain, while the corresponding eigenvalue at K does not show a big variation. This behavior can be explained by analyzing the character of the KS wavefunction at K v and at G v (see bottom panel of Figure 3). While the main contribution of K v is due to p z N orbitals, G v has a large contribution from p x;y (VBM and VBM-1 are degenerated) orbitals of N atoms. The analysis of the corresponding states in AlN (see Figure 3, top panel) reveals similar contributions for K v and G v , but the more localized nature of the orbitals implies a smaller sensitivity of the bonds to the applied strain and, in particular, a minor energetic variation of G v in AlN under compressive strain, with respect to the case of GaN. We report in Figure 4, for 2D AlN and GaN, the dependence of valence band maximum (VBM), conduction band minimum (CBM), vacuum potential (left panel); IP, EA, WF (center panel), and electronic gaps (right panel). In both materials, the vacuum potential, as well as CBM and VBM, decreases (increases) when increasing tensile (compressive) strain (left panel). Looking at the central panel, we note that the IP of GaN increases (decreases) with increasing tensile (compressive) strain, showing similar behavior to what observed for the WF of X-enes. The two different slopes between À10% to À4% and between À4% to 10% are related to the change of VBM from G v to K v , as discussed before.
It is worth to notice that the corresponding IP of AlN shows an opposite trend, increasing (decreasing) when the compressive (tensile) strain increases. This can be explained, again, in terms of the stronger localization of  the electronic charge near the nitrogen atoms that occurs in AlN, with respect to the case of GaN, due to the larger difference in the electronegativity between anion and cation. The EA is given by the position of the conduction bands bottom, and both materials show a similar behavior: EA increases (decreases) for increasing tensile (compressive) strain. The CBM is always at G and from the analysis of the KS wave functions, it is due to a hybrid of s states of the cation (Al/Ga) and anion (N), where the former gives a larger contribution. Also, for GaN, a small contribution from the d states appears. This analysis explains why the slope of CBM and then of EA is steeper in GaN with respect to AlN. The dependence of WF and gaps from strain is then a direct consequence of the IP and EA behavior.

TMD: MoTe 2 and MoS 2
Upon the application of strain, the electronic properties of TMD show a more complicate behavior with respect to those of X-enes and nitrides. Moreover, despite the similarity in their atomic structures, the electronic properties of MoS 2 and MoTe 2 monolayers and bilayers have similar but not completely equivalent trends. 64,65 MoTe 2 and MoS 2 monolayer. In agreement with the existing literature, the two MoX 2 -ML studied here exhibits, at zero strain, a direct bandgap with VBM and CBM located at the six K points corners of the BZ (see the central panels of   calculations performed at the same level of theoretical approximation. 66 Beyond the K points also G and Q (where a minimum of the CB is visible in the band structure, along the G À K direction) critical points have to be mentioned, because they are energetically close to K v and K c and, under strain or in multilayers, Q c may become the global CBM and G v the global VBM. For all the group-VI TMD, it is well-known that the edges of the conduction and valence bands are composed predominantly by d metal (M) and p chalcogen (X) orbitals. 67 However, the types (symmetries) and contributions of these orbitals vary with the chosen kpoint and with the constituent atoms. As shown in Figure 6 at K, the VB (CB) is mainly described by d xy ; d x 2 Ày 2 (d z 2 ) M-orbitals with a small contribution of p x;y X-orbitals. At G, the VB has mainly a M-d z 2 character with a small X-p z contribution.
The analysis of Q c shows a similar composition of K v with minor contribution from d z 2 and p x;y;z orbitals. 68 Uniform tensile strain causes a general reduction of the direct gap (at K), which can be understood, using a simple tightbinding picture, as due to a minor overlap among orbitals. Nevertheless, because of the different orbital composition, the application of strain acts in a different way at the various high-symmetry k-points. In particular, due to the dominant d xy ; d x 2 Ày 2 contribution, a downward shift of K v is observed, while G v , for its d z 2 main composition, remains essentially unshifted. 69 This effect induces a direct-toindirect bandgap crossover (G v À K c ) when tensile strain is applied, which occurs at e < 1ð*5Þ% for MoS 2 (MoTe 2 ) monolayer. Indeed, the extreme delocalized nature (see, for instance, Figure 7) of the p-orbitals of Te with respect to those of sulfur makes MoTe 2 less sensitive to strain. This means that smaller amounts of tensile strain are needed to shift down K v in MoS 2 , as evident in the right panels of Figure 5.
Interestingly, the behavior of MoS 2 and MoTe 2 under uniform compressive strain is very different. In MoS 2 -ML, the CBM moves from K to Q for strain above À2%, hence inducing the formation of an indirect gap K v À Q c . On the other hand, in MoTe 2 -ML, when compressive strain increases, the two CB minima, along the L and S directions decrease their energies, while the valence band at M increases in energy. These two facts lead to the first transition from a direct gap K c À K v to an indirect gap Q c À K v at about À2% of strain and to a second transition to another indirect gap, Q c À M v , for strain of about À5%. From the analysis of the wavefunctions at M v , it results that contributions come from d x 2 Ày 2 ,d z 2 orbitals of the metal atoms and more dominant contributions come from p x;y;z orbitals of chalcogen atoms. Again, due to the more delocalized nature of the p orbitals of Te with respect to S, M v is much more sensitive to compressive strain in MoTe 2 than in MoS 2 (left panel of Figure 8). The calculated WF, IP, EA, and electronic gap trends with strain for TMD monolayers are shown in Figure 8. VBM, CBM, and vacuum potential curves (left panel) are qualitatively similar to those of 2D nitrides, but we note that the slopes are very different. In particular, the vacuum potential decreases more rapidly in MoTe 2 with respect to the other cases. For MoS 2 -ML, the IP (central panel in Figure 8) decreases when compressive strain is applied with an almost linear behavior due to the fact that the VBM remains K v . Also, in MoTe 2 , the IP decreases by applying a compressive strain, but at approximately À5%, there is a change in the slope, caused by the change of the VBM from K v to M v .  For what concerns tensile strain, we see that the IP of MoS 2 decreases by increasing tensile strain. This can be explained in the following way: from very small strain values (less than 1%), the VBM moves from K v to G v and it remains almost unaffected by the strain (see the blue dots in the left panel of Figure 8); in this way, we can attribute the IP variation to the decreasing behavior of the vacuum potential. For MoTe 2 (black triangles, center panel in Figure 8), IP first increases until it reaches a maximum around þ5% and then begins to decrease due to the fact that the VBM shift from K to G at higher tensile strains.
For both TMD-ML, EA decreases when reducing the lattice parameter up to the crossover of the CBM from K c to Q c . Then, EA begins to increase with a very small rate when Q c becomes the global conduction band minimum (central panel of Figure 8).
From the analysis of the data shown in the central panel of Figure 8, and similar to the case of nitrides, it is clear that the behavior of the WF depends on both IP and EA curves and shows an intermediate behavior. We can say that the WF and EA values for the two TMD-ML have a qualitatively similar behavior: an increase with a change in the slope that occurs when the direct to indirect bandgap crossover happens for low compressive strain. The linearity of EA is due to the smoothness of the vacuum potential and   much more symmetrical behavior of MoTe 2 for tensile and compressive strain.
Beyond the monolayer. TMD bilayers and multilayers are known to be indirect bandgap semiconductors. MoS 2 (MoTe 2 ) shows an indirect gap G v À K c of 1.1 eV (K v À Q c of 0:93 eV) in partial agreement with previous results 47 (see the central panels of Figure 9).
The p X-orbital composition at Q c and G v plays an important role in the crossover from the direct to indirect bandgap from monolayer to bulk. Indeed, the close distance between X-p orbitals from neighboring layers leads to large hopping, which changes the energy of Q c and G v substantially. 67 For this reason, the bandgap is very sensitive to small changes in the interlayer distance, thus dependent on the approximated vdW functional used.
In this regard, it is worth to mention that for MoTe 2 bulk, we specifically tested how our simulated structural parameters (a ¼ 3.54Å, Mo-Mo vertical distance between two different layers 7.33Å) are obtained using the selected vdW functional, reproduce the experimental data (a ¼ 3.54 A, Mo-Mo distance between two different layers 7.33Å) reported by Knop and MacDonald,70 finding a very good agreement.
Moreover, always for MoTe 2 -BL, we also tested other two vdW functionals taken from the study by Grimme 71 and Barone et al., 72 respectively. In both cases, the final optimized structures are very similar to the one reported in Table 1, with lattice parameter a and Te-Te intralayer vertical distance, which changes less than 0.3%. Only using the second functional, we observe a reduction of 2.6% of the interlayer vertical distance. The electronic band structure does not show relevant changes near the gap region with the largest downshift of the VBM at G of 0.05 eV.
In Figure 7, we show the charge density plot of the Q c state for MoS 2 (left) and MoTe 2 (right) bilayers, where the important contribution coming from p-orbitals of chalcogen atoms is evident. Similar to what discussed above, the very delocalized nature of these orbitals for the case of Te is evident from the figure. MoS 2 and MoTe 2 bilayers show, already at zero strain, different values of VBM and CBM due to the different hopping between p of the two adjacent layers, as shown in Figure  7. When tensile strain is applied in both BL, K c and K v go down in energy in agreement with the monolayer behavior (right panels of Figure 9). For compressive strain, similar to the monolayers, Q c goes down in energy and, for MoTe 2 , M v becomes the VBM when À4% strain is applied (left panel).
In other words, trends previously discussed for monolayers remain qualitative similar for bilayers (see Figure 9), but different values of strain are needed to obtain changes in VBM and CBM with respect to the MLs. Moreover, MoTe 2 -BL results are more sensitive to strain with respect to MoS 2 -BL. This is also highlighted in Figure10, where the IP curves of the monolayer and bilayer of MoS 2 and MoTe 2 are reported, respectively, in the top and bottom panels.

Role of many-body effects
Despite the DFT electronic properties of 2D materials that are affected by the well-known bandgap problem, previous works have suggested that their behavior under strain reproduces, at least qualitatively, the one obtained with more refined excited state methods, such as GW . 65 In particular, for MoTe 2 , it has been shown 73 that the response to strain of CBM and VBM studied within DFT plus hybrid functional HSE06 qualitatively agree. To confirm these results, we have calculated, for MoTe 2 -ML, how CBM and VBM change with strain, using the perturbative G 0 W 0 method. The obtained G 0 W 0 bandgap is 1.73 eV for the unstrained structure. This value is in good agreement with that one reported by Robert et al. 74 and Rasmussen and Thygesen 75 of 1.72 eV.
Because of the heaviness of GW calculations, we have limited our analysis to zero strain and to one value of compressive (À2%) and tensile (2%) strain. In Figure 11, we report the energy values of the CBM (red) and VBM (blue) at both DFT and G 0 W 0 levels of approximation. The VBM and the CBM are at K for 0 and þ2% tensile strain, but the CBM moves at Q for À2% strain. While the KS eigenvalues (dashed lines) for K c and K v follow a trend with strain qualitatively similar to the correct quasiparticle one (solid lines), some deviation occurs for Q c . This is clearly due to the fact that the quasiparticle correction of the state at Q c is quite different from that of the state at K c , due to the different contributions of atomic orbitals composing them. Our analysis of many-body effects is clearly limited, having considered only three values of strain, but it is clear that some caution has to be taken when trends with strain of the electronic gaps and work functions are extracted at the DFT level if a change in the character of the valence and conduction states near the gap occurs.

Conclusions
We have presented a systematic study of the effect of uniform biaxial compressive and tensile strain on the electronic properties of several 2D materials ranging from semimetallic ones, like the X-enes, to semiconductors like nitrides and TMD.
Each material shows its own behavior strictly linked to the changes in the band structures and to the orbital character of the CBM and VBM. In all monolayers, the WF always increases with the lattice parameter, while the IP and EA show material-dependent trends.
The slope of the IP versus strain is usually positive, but exceptions arise for AlN (very small negative slope) and MoS 2 (where the right and left derivatives are different). The bilayers appear to be much more sensitive to strain than the isolated layers, hence, the curves are similar only in some regions. The difference in behaviors for nitrides and TMD is explained in terms of the localization of the VBM and CBM wave functions. We have also shown that the choice of the XC potential while affecting the absolute value of the WF does not change the general behavior of the electronic properties with respect to the strain. Finally, we have demonstrated that special care has to be used when the character of the VBM or CBM changes with the strain, as the DFT and G 0 W 0 trends may not be equivalent.