Structure–Property Relationships in Transition Metal Dichalcogenide Bilayers under Biaxial Strains

This paper reports a Density Functional Theory (DFT) investigation of the electron density and optoelectronic properties of two-dimensional (2D) MX2 (M = Mo, W and X = S, Se, Te) subjected to biaxial strains. Upon strains ranging from −4% (compressive strain) to +4% (tensile strain), MX2 bilayers keep the same bandgap type but undergo a non-symmetrical evolution of bandgap energies and corresponding effective masses of charge carriers (m*). Despite a consistency regarding the electronic properties of Mo- and WX2 for a given X, the strain-induced bandgap shrinkage and m* lowering are strong enough to alter the strain-free sequence MTe2, MSe2, MS2, thus tailoring the photovoltaic properties, which are found to be direction dependent. Based on the quantum theory of atoms in molecules, the bond degree (BD) at the bond critical points was determined. Under strain, the X-X BD decreases linearly as X atomic number increases. However, the kinetic energy per electron G/ρ at the bond critical point is independent of strains with the lowest values for X = Te, which can be related to the highest polarizability evidenced from the dielectric properties. A cubic relationship between the absolute BD summation of M-X and X-X bonds and the static relative permittivity was observed. The dominant position of X-X bond participating in this cubic relationship in the absence of strain was substantially reinforced in the presence of strain, yielding the leading role of the X-X bond instead of the M-X one in the photovoltaic response of 2D MX2 material.

In a previous work [28], we investigated the layered-dependent structural, electronic, and optical properties of MX 2 homo-and heterostructures by DFT calculations. The quantum theory of atoms in molecules [29] was used to process the electron density in order to correlate electronic interactions and macroscopic optical properties. We found that the static relative permittivity and the weighted bond degree summation are linked by a cubic relation and that the layered-dependent electronic and optical properties are mainly attributed to the interlayer X-X bonds. Furthermore, it has been reported in literature that the application of strains allows one to finely tune the transition metals dichalcogenides (TMD) electronic properties [30][31][32][33][34][35][36]. For instance, while uniaxially strained, the MoS 2 monolayer changes from direct to indirect bandgap, which opens up the possibility of designing TMD through strain engineering. In addition, whereas conventional materials hardly bear strains exceeding a few percent, MoS 2 has been shown to be able to withstand strains above 11% [37].
A great number of first principle calculations have already been reported on the effect of strain on the electronic properties of 2D semiconducting TMD [36,[38][39][40][41][42]. However, these works almost all focus on monolayers. As for the effect of strain on the optical properties, it has only been investigated for monolayers. In a comparative investigation on MoS 2 monolayers, Carrascoso et al. [43] showed that uniaxial strains have a weaker effect on the materials properties than biaxial ones. This type of comparative study has never been carried out on bilayer TMD, but we expect a similar trend. Hence, building upon our previous work, we investigated for the first time the effect of biaxial strains on both electronic and optical properties and the contribution of the interlayer van der Waals interactions in the optical properties in the bilayered MX 2 compounds.

Materials and Methods
DFT [44,45] calculations were carried out by a full-potential linear augmented plane wave method (FP-LAPW) as implemented in the program WIEN2k [46]. The generalized gradient approximation level of theory was applied with the Wu-Cohen (WC) functional [47]. During the optimization of the atomic position, the convergence criteria were set to 10 −5 Ry and 1 mRy/Bohr for the energy and forces, respectively. The R mt K max parameter was set to 7. Besides, the first Brillouin zone was sampled with 1500 k-points that were mesh selected according to the Monkhorst-Pack algorithm [48]. Although hybrid range-separated functionals are now recognized as a standard for obtaining an accurate description of chemical systems, the choice of the WC functional was made after systematic tests and comparisons with available data. It appears that this functional gives very reliable structural parameters compared to experimental ones. Regarding band structures, this functional also yields decent results, which may be attributed to the fact that we are investigating chalcogenides, for which the energy gap is rarely wrongly zeroed. In addition, considering the huge amount of calculations achieved in this work, we concluded that using a more elaborated functional, such as a hybrid one, would have led to such a computational cost that it could have jeopardized the achievement of our objectives.
The investigated MX 2 (M = Mo, W; X = S, Se, Te) bilayers consisting of 4 × 4 × 1 supercells and the density of states near the Fermi level (where MoS 2 is taken as an example) are depicted in Figure 1. A 20 Å vacuum thickness was added atop to separate free surfaces, hence avoiding interaction between periodic images. The structures were then relaxed. Subsequently, both compressive and tensile in-plane biaxial strains were applied with values ranging from -4% to +4% by steps of 2%. Negative deformations stand for compressive strains, whereas positive ones stand for tensile strains. After structure relaxation, electronic band structures and optical properties (relative permittivity, absorption coefficient, extinction coefficient, and refractive index) were calculated for each of the six MX 2 bilayers. The application of a biaxial strain on chalcogenide bilayers aims at reproducing the epitaxial strain on the absorber layer in real devices. Nonetheless, our model is limited by the absence of a substrate, which does not allow us to investigate the band alignments in the device. Hence, the energy discontinuities at the band edges that serve as the basis for controlling transport properties were not characterized. Further, the electron densities obtained from WIEN2k calculations were processed with the Critic2 package [49], which implements the quantum theory of atoms in molecules [29,50], from which the total, kinetic, and potential energy densities at the bond critical points (BCPs) were obtained. The bond degree at each BCP was then calculated from the total energy density and the electron density [51].

Electronic and Optical Properties
The band structures of the six investigated MX2 bilayers under biaxial strains (from compressive to tensile ones) are depicted in Figure 2. For the unstrained bilayers (null strain), except for WTe2, all the bilayers bear an indirect bandgap that decreases from S to Se and then to Te. The fact of applying strains, either compressive or tensile, does not change the type of bandgap (again, except in the case of WTe2, which becomes indirect), but the shape of the bands can be substantially modified, potentially implying a change in the valence band minimum (VBM) and conduction band maximum (CBM); the k-points implied in the valence-to-conduction transitions are hence also changed accordingly. The bandgap evolution with respect to the applied strains is depicted in Figure 3a. It shows that, except for Te, the bandgaps in the WX2 bilayers are larger than in the MoX2 ones for a given chalcogen atom. Moreover, increasing the tensile strain leads to a bandgap shrinkage, which increases along the MTe2, MSe2, and MS2 sequence, irrespective of the metal atom. These tensile strain effects are similar to those observed when the materials go from a bulk state to a monolayer one: a blue shift of the energy bandgap is observed, which is attributed to quantum confinement [52]. Except for Te, for a given chalcogen atom, the shrinkage decrease is higher when the W metal is concerned. By contrast, the bandgap does not change in a systematic way as the compressive strain increases but clearly depends on the nature of both the metal atom and the chalcogen one. For WSe2 and MoTe2, the bandgap decreases when the compressive strain increases, whereas it increases for MoS2 and WS2. For the remaining compounds (MoSe2 and WTe2), the bandgap first increases and then decreases when the compressive strain increases. These results agree with those reported in the literature on bilayer MX2 under both uniaxial and biaxial strains [35,53]. As reported by Carrascoso et al. for monolayers, the uniaxial strains have a weaker effect on the bilayer material properties than the biaxial ones [43]. Further, the electron densities obtained from WIEN2k calculations were processed with the Critic2 package [49], which implements the quantum theory of atoms in molecules [29,50], from which the total, kinetic, and potential energy densities at the bond critical points (BCPs) were obtained. The bond degree at each BCP was then calculated from the total energy density and the electron density [51].

Electronic and Optical Properties
The band structures of the six investigated MX 2 bilayers under biaxial strains (from compressive to tensile ones) are depicted in Figure 2. For the unstrained bilayers (null strain), except for WTe 2 , all the bilayers bear an indirect bandgap that decreases from S to Se and then to Te. The fact of applying strains, either compressive or tensile, does not change the type of bandgap (again, except in the case of WTe 2 , which becomes indirect), but the shape of the bands can be substantially modified, potentially implying a change in the valence band minimum (VBM) and conduction band maximum (CBM); the k-points implied in the valence-to-conduction transitions are hence also changed accordingly. The bandgap evolution with respect to the applied strains is depicted in Figure 3a. It shows that, except for Te, the bandgaps in the WX 2 bilayers are larger than in the MoX 2 ones for a given chalcogen atom. Moreover, increasing the tensile strain leads to a bandgap shrinkage, which increases along the MTe 2 , MSe 2 , and MS 2 sequence, irrespective of the metal atom. These tensile strain effects are similar to those observed when the materials go from a bulk state to a monolayer one: a blue shift of the energy bandgap is observed, which is attributed to quantum confinement [52]. Except for Te, for a given chalcogen atom, the shrinkage decrease is higher when the W metal is concerned. By contrast, the bandgap does not change in a systematic way as the compressive strain increases but clearly depends on the nature of both the metal atom and the chalcogen one. For WSe 2 and MoTe 2 , the bandgap decreases when the compressive strain increases, whereas it increases for MoS 2 and WS 2 . For the remaining compounds (MoSe 2 and WTe 2 ), the bandgap first increases and then decreases when the compressive strain increases. These results agree with those reported in the literature on bilayer MX 2 under both uniaxial and biaxial strains [35,53]. As reported by Carrascoso et al. for monolayers, the uniaxial strains have a weaker effect on the bilayer material properties than the biaxial ones [43].  All the aforementioned effects undoubtedly have an impact on the valence-to-conduction electron transitions upon irradiation and on the excited electrons' mobility due to band curvature changes. Indeed, the best voltage and photocurrent amplitudes are dictated by the proper balancing between the bandgap and absorbed photon energy [54]. Figure 3b plots the corresponding electron and hole effective masses (me* and mh*, respectively) at the respective conduction and valence band edges. For each fixed X atom, MoX2 has a comparatively higher m* than WX2. Specifically, the me* and mh* along the Γ-K direction become smaller as the tensile strain scales up, except for MoTe2 and WTe2 in the range of [0%; +2%] for which they become larger. However, in the range of [0%; −2%], the me* and mh* along either Γ-Λ or K-Λ decrease and increase, respectively. As the compressive strain intensifies, the me* and mh* along the K-Λ direction decrease again. It is worth mentioning that the me* and mh* of MoTe2 do not follow this pattern due to the change of the electron transition path from the K-Λ to the M-Λ direction. All the above information confirms an influence of the strain on the band edges and curvatures, which can be deep enough to shift the sequence of bandgap energies and m* in MX2 for different X, thus providing possibilities of customizing the photovoltaic properties. Hence, considering the ideal bandgap and favorable m*, an advantageous electron excitation and transportation can be anticipated when a certain range of compressive strain is applied.
In the following, the optical properties were analyzed in the in-plane (xx) and out-ofplane (zz) directions. The calculated absorption coefficients and refractive indexes of the  All the aforementioned effects undoubtedly have an impact on the valence-to-co duction electron transitions upon irradiation and on the excited electrons' mobility due band curvature changes. Indeed, the best voltage and photocurrent amplitudes are d tated by the proper balancing between the bandgap and absorbed photon energy [5 Figure 3b plots the corresponding electron and hole effective masses (me* and mh*, respe tively) at the respective conduction and valence band edges. For each fixed X atom, Mo has a comparatively higher m* than WX2. Specifically, the me* and mh* along the Γ-K d rection become smaller as the tensile strain scales up, except for MoTe2 and WTe2 in t range of [0%; +2%] for which they become larger. However, in the range of [0%; −2%], t me* and mh* along either Γ-Λ or K-Λ decrease and increase, respectively. As the compre sive strain intensifies, the me* and mh* along the K-Λ direction decrease again. It is wor mentioning that the me* and mh* of MoTe2 do not follow this pattern due to the change the electron transition path from the K-Λ to the M-Λ direction. All the above informati All the aforementioned effects undoubtedly have an impact on the valence-to-conduction electron transitions upon irradiation and on the excited electrons' mobility due to band curvature changes. Indeed, the best voltage and photocurrent amplitudes are dictated by the proper balancing between the bandgap and absorbed photon energy [54]. Figure 3b plots the corresponding electron and hole effective masses (m e * and m h *, respectively) at the respective conduction and valence band edges. For each fixed X atom, MoX 2 has a comparatively higher m* than WX 2 . Specifically, the m e * and m h * along the Γ-K direction become smaller as the tensile strain scales up, except for MoTe 2 and WTe 2 in the range of [0%; +2%] for which they become larger. However, in the range of [0%; −2%], the m e * and m h * along either Γ-Λ or K-Λ decrease and increase, respectively. As the compressive strain intensifies, the m e * and m h * along the K-Λ direction decrease again. It is worth mentioning that the m e * and m h * of MoTe 2 do not follow this pattern due to the change of the electron transition path from the K-Λ to the M-Λ direction. All the above information confirms an influence of the strain on the band edges and curvatures, which can be deep enough to shift the sequence of bandgap energies and m* in MX 2 for different X, thus providing possibilities of customizing the photovoltaic properties. Hence, considering the ideal bandgap and favorable m*, an advantageous electron excitation and transportation can be anticipated when a certain range of compressive strain is applied.
In the following, the optical properties were analyzed in the in-plane (xx) and out-ofplane (zz) directions. The calculated absorption coefficients and refractive indexes of the investigated MX 2 compounds are shown in Figures 4 and 5, respectively. Irrespective of the compound and strain, the absorption threshold and refraction peak occur at a lower energy in the xx direction than in the zz one, with a slight, gradual shift towards higher energies from +4% to −4% strain in the case of the xx-direction; in the zz-direction, the curves are indistinguishable. Thus, for a given compound, a tensile strain allows for lowering the absorption threshold. This result agrees with the decreasing bandgap observed under tensile strain and is similar to that observed by a photoluminescence spectroscopy experiment, as the materials' size decreases from the bulk to monolayer [3,4,55]. Irrespective of the metal atom and strain, the refraction peak intensity increases, and the absorption edge value decreases as the chalcogen atomic number increases. investigated MX2 compounds are shown in Figures 4 and 5, respectively. Irrespective of the compound and strain, the absorption threshold and refraction peak occur at a lower energy in the xx direction than in the zz one, with a slight, gradual shift towards higher energies from +4% to −4% strain in the case of the xx-direction; in the zz-direction, the curves are indistinguishable. Thus, for a given compound, a tensile strain allows for lowering the absorption threshold. This result agrees with the decreasing bandgap observed under tensile strain and is similar to that observed by a photoluminescence spectroscopy experiment, as the materials' size decreases from the bulk to monolayer [3,4,55]. Irrespective of the metal atom and strain, the refraction peak intensity increases, and the absorption edge value decreases as the chalcogen atomic number increases.   The relative permittivity function ε (ω) is strongly related to the band structure and characterizes collective excitations close to the Fermi level [56]. The calculated dielectric functions for the six MX2 bilayers, subjected or not to strains, are depicted in Figure 6, with the real part ε1 (ω) being shown on the left panel and imaginary one ε2 (ω) on the right panel. For all the compounds, ε1 (ω) becomes negative above around 5-6 eV, which means that the compounds exhibit a metallic behavior above the photon energy thresholds [57]. At the frequency limit ω = 0, ε1 (0) corresponds to the static relative permittivity. For all the compounds, higher corresponding values were obtained in the xx-direction than in the zz-direction, and they increase from compressive to tensile strain, whereas they decrease in the zz-direction. Irrespective of the metal atom, ε1 (0) increases with the chalcogen atomic number. Hence, based on the Penn model [58], which defines ε1 (0) as 1 (0) ≈ 1 + ( ℏ ) 2 , the bandgap should decrease from S to Te, which is indeed observed in Figure   3a. The highest static relative permittivity value was obtained for MoTe2, indicating a higher polarizability for this bilayer. These values are further improved in the xx-and zzdirection when the compound undergoes tensile and compressive strains, respectively. For all the compounds, when compared with ε2 (ω) in the xx-direction, we observe that ε2 (ω) in the zz-direction tends to decrease and that its peak maximum shifts towards higher incident photon energies. These results indicate a decrease in the ability of the compounds to absorb light in this direction. The relative permittivity function ε (ω) is strongly related to the band structure and characterizes collective excitations close to the Fermi level [56]. The calculated dielectric functions for the six MX 2 bilayers, subjected or not to strains, are depicted in Figure 6, with the real part ε 1 (ω) being shown on the left panel and imaginary one ε 2 (ω) on the right panel. For all the compounds, ε 1 (ω) becomes negative above around 5-6 eV, which means that the compounds exhibit a metallic behavior above the photon energy thresholds [57]. At the frequency limit ω = 0, ε 1 (0) corresponds to the static relative permittivity. For all the compounds, higher corresponding values were obtained in the xx-direction than in the zz-direction, and they increase from compressive to tensile strain, whereas they decrease in the zz-direction. Irrespective of the metal atom, ε 1 (0) increases with the chalcogen atomic number. Hence, based on the Penn model [58], which defines ε 1 (0) as 1 (0) ≈ 1 + ω E g 2 , the bandgap should decrease from S to Te, which is indeed observed in Figure 3a. The highest static relative permittivity value was obtained for MoTe 2 , indicating a higher polarizability for this bilayer. These values are further improved in the xxand zz-direction when the compound undergoes tensile and compressive strains, respectively. For all the compounds, when compared with ε 2 (ω) in the xx-direction, we observe that ε 2 (ω) in the zz-direction tends to decrease and that its peak maximum shifts towards higher incident photon energies. These results indicate a decrease in the ability of the compounds to absorb light in this direction.

Electron Density Analysis
In the realm of the quantum theory of atoms in molecules (QTAIM) [29,50] the key fields are the electron density, and especially its Laplacian, from which numerous parameters were derived, which enables the characterization of the bonding between atoms in molecules and crystals [59]. Among these parameters, the bond degree (BD) at the bond critical point BD = H b /ρ b , where H is the total energy density, i.e., the sum of the potential V and kinetic G energy densities, and ρ is the electron density at the bond critical point (b), measures the degree of covalence (BD < 0) or softening (BD > 0) of the interatomic bonding [51]. In other words, covalent bonds are characterized by large, negative values of BD, whereas closed-shell interactions (ionic and van der Waals interactions) are characterized by positive BD values. Figures 7 and 8 depict the bond degree of both the M-X and X-X bonds for the bilayers subjected to strains. In addition, the evolution of the M-X and X-X bond lengths under strains are also depicted in Figure 7. Unsurprisingly, for both the M-X and X-X bonds, the BL increases with the chalcogen atomic number, irrespective of the metal atom. For a given chalcogen atom, the Mo-X and W-X bond lengths are the same for each applied strain, which can be explained by the similar value of the Mo and W covalent radii (145 pm and 146 pm, respectively). The M-X bond lengths linearly increase when the strain varies from −4% to +4%. In the case of the X-X bond lengths, a slight difference is noticeable depending on whether Mo or W is bonded to the chalcogen atom. This slight difference is also observed in the corresponding bond degrees, which, contrary to the bond lengths, decrease when the chalcogen atomic number increases. This decrease reflects the lowering of the van der Waals character of the interatomic interaction. Regarding the evolution with the strain, both the X-X bond lengths and bond degrees are nearly constant. For the M-X bonds, the BDs increase from compressive to tensile strains, the absolute values of which are increasingly large for S, Te, and Se.

Electron Density Analysis
In the realm of the quantum theory of atoms in molecules (QTAIM) [29,50] fields are the electron density, and especially its Laplacian, from which numerous eters were derived, which enables the characterization of the bonding between a molecules and crystals [59]. Among these parameters, the bond degree (BD) at t critical point BD = Hb/ρb, where H is the total energy density, i.e., the sum of the p V and kinetic G energy densities, and ρ is the electron density at the bond critic (b), measures the degree of covalence (BD < 0) or softening (BD > 0) of the inte bonding [51]. In other words, covalent bonds are characterized by large, negativ of BD, whereas closed-shell interactions (ionic and van der Waals interactions) a acterized by positive BD values. Figures 7 and 8 depict the bond degree of both the M-X and X-X bonds for the subjected to strains. In addition, the evolution of the M-X and X-X bond length strains are also depicted in Figure 7. Unsurprisingly, for both the M-X and X-X bo BL increases with the chalcogen atomic number, irrespective of the metal atom given chalcogen atom, the Mo-X and W-X bond lengths are the same for each strain, which can be explained by the similar value of the Mo and W covalent ra pm and 146 pm, respectively). The M-X bond lengths linearly increase when th varies from −4% to +4%. In the case of the X-X bond lengths, a slight difference is ble depending on whether Mo or W is bonded to the chalcogen atom. This slight di is also observed in the corresponding bond degrees, which, contrary to the bond decrease when the chalcogen atomic number increases. This decrease reflects the l of the van der Waals character of the interatomic interaction. Regarding the evolut the strain, both the X-X bond lengths and bond degrees are nearly constant. For bonds, the BDs increase from compressive to tensile strains, the absolute values o are increasingly large for S, Te, and Se.   According to Figure 8, there is no clear relation between BD and |V|/G values of M-X bonds. By contrast, a linear relationship can be evidenced for the X-X ones, the slope being the same for all the compounds (see Figure 8g) with or without strain. However, irrespective of the bond and the compound, the evolution of both BD and |V|/G values with respect to strain seems to be weak. In order to better evaluate this evolution, the G/ρ value, which corresponds to the slope of the line passing through the point of interest and that of the coordinates |V|/G = 1 and BD = 0 (see [60]), was determined for each bond type in each bilayer with and without strains. The results are gathered in Table 1. One can see the following: (i) the G/ρ values for each bond are almost independent of strain; (ii) the G/ρ values of both X-X and M-X bonds mainly depend on the X atom, the influence of the M one being very weak; (iii) the G/ρ values decrease when the atomic number of the X atom implied in the bond increases; (iv) for X= Se and Te, the G/ρ values of M-X and X-X bonds are very close to each other, and for X = S, the G/ρ values of M-X bonds are higher than those of the X-X ones. These results agree with the inverse relation between the kinetic energies per electron G/ρ at the bond critical point and the bond polarizability, as proposed by Yang et al. [60]. Indeed, the larger chemical softness of Te compared to that of Se and S, and the larger one of Mo compared to that of W, should correspond to larger polarizability. This can be related to the highest polarizability evidenced in the previous section for MoTe2 from static relative permittivity values. According to Figure 8, there is no clear relation between BD and |V|/G values of M-X bonds. By contrast, a linear relationship can be evidenced for the X-X ones, the slope being the same for all the compounds (see Figure 8g) with or without strain. However, irrespective of the bond and the compound, the evolution of both BD and |V|/G values with respect to strain seems to be weak. In order to better evaluate this evolution, the G/ρ value, which corresponds to the slope of the line passing through the point of interest and that of the coordinates |V|/G = 1 and BD = 0 (see [60]), was determined for each bond type in each bilayer with and without strains. The results are gathered in Table 1. One can see the following: (i) the G/ρ values for each bond are almost independent of strain; (ii) the G/ρ values of both X-X and M-X bonds mainly depend on the X atom, the influence of the M one being very weak; (iii) the G/ρ values decrease when the atomic number of the X atom implied in the bond increases; (iv) for X= Se and Te, the G/ρ values of M-X and X-X bonds are very close to each other, and for X = S, the G/ρ values of M-X bonds are higher than those of the X-X ones. These results agree with the inverse relation between the kinetic energies per electron G/ρ at the bond critical point and the bond polarizability, as proposed by Yang et al. [60]. Indeed, the larger chemical softness of Te compared to that of Se and S, and the larger one of Mo compared to that of W, should correspond to larger polarizability. This can be related to the highest polarizability evidenced in the previous section for MoTe 2 from static relative permittivity values. According to Gatti [59], an atomic expectation value results from the sum of bond contributions. As we did in previous works [28,61], a relationship was searched for between the bond degrees summation and the relative permittivity under zero frequency ε 1 (0) along the zz-direction. As two types of bonds coexist in the structures, namely M-X and X-X, the summation can be written as h|BD| M-X +k|BD| X-X , with h and k the parameters to be fitted. Fitting this expression via the equation, the maximum coefficient of determination R 2 is obtained by adjusting h and k at each equation order. Irrespective of the strain, the most accurate description of the relationship between bond degree and static relative permittivity is given by a cubic equation (see Figure 9a), namely, n = 3. The absolute BD summation and the ε 1 (0) are inversely related. The best fit at n = 3, for which R 2 = 0.945, was obtained for h/k = 0.05, in comparison to the best fit at n = 1 and n = 2, for which R 2 = 0.849 and R 2 = 0.881, respectively, both obtained for h/k = 0.00. These results indicate that in strain-modified MX 2 bilayers, the X-X bonds are overwhelmingly contributing to the dielectric properties. By contrast, the fitting result under no strain, as seen in Figure 9b, shows a perfect cubic relationship between the bond degree summation and static relative permittivity with the R 2 = 0.997 for h/k = 0.3. This is coherent with our previous observation in the absence of strains [28], where both types of bonds were found to participate in achieving this cubic relationship, although the X-X bonds were found to contribute more than the M-X ones. The profound decrease of h/k ratio of MX 2 from no applied strain to added strain highlights the pronouncing role of the X-X bonds in responding to the external photoelectric field, as proven by the manifest dependence of the absorption and refractive properties on the strain as X varies.
According to Gatti [59], an atomic expectation value results from the sum of bond contributions. As we did in previous works [28,61], a relationship was searched for between the bond degrees summation and the relative permittivity under zero frequency ε1 (0) along the zz-direction. As two types of bonds coexist in the structures, namely M-X and X-X, the summation can be written as h|BD|M-X+k|BD|X-X, with h and k the parameters to be fitted. Fitting this expression via the equation, the maximum coefficient of determination R 2 is obtained by adjusting h and k at each equation order. Irrespective of the strain, the most accurate description of the relationship between bond degree and static relative permittivity is given by a cubic equation (see Figure 9a), namely, n = 3. The absolute BD summation and the ε1 (0) are inversely related. The best fit at n = 3, for which R 2 = 0.945, was obtained for h/k = 0.05, in comparison to the best fit at n = 1 and n = 2, for which R 2 = 0.849 and R 2 = 0.881, respectively, both obtained for h/k = 0.00. These results indicate that in strain-modified MX2 bilayers, the X-X bonds are overwhelmingly contributing to the dielectric properties. By contrast, the fitting result under no strain, as seen in Figure  9b, shows a perfect cubic relationship between the bond degree summation and static relative permittivity with the R 2 = 0.997 for h/k = 0.3. This is coherent with our previous observation in the absence of strains [28], where both types of bonds were found to participate in achieving this cubic relationship, although the X-X bonds were found to contribute more than the M-X ones. The profound decrease of h/k ratio of MX2 from no applied strain to added strain highlights the pronouncing role of the X-X bonds in responding to the external photoelectric field, as proven by the manifest dependence of the absorption and refractive properties on the strain as X varies. Inset in (a) depicts the first-, second-and third-order fittings of R 2 vs. h/k. Inset in (b) depicts the third-order fitting of R 2 vs. h/k ratio. Arrows in between symbols "−" and "+" in (a) are used to represent the evolution from the biaxial compressive (−4%, −2%) strains to the tensile (+2%, +4%) ones.

Conclusions
The influence of biaxial strain on optoelectronic properties and electron density of transition metal dichalcogenides MX2 (M = Mo, W and X = S, Se, Te) bilayers were thoroughly investigated for the first time using DFT calculations. When subjected to a strain going from a compressive (−4%, −2%) to tensile (+2%, +4%) one, the 2D materials' band structures and their corresponding effective masses of charge carriers undergo non-symmetrical changes as compression and tension are concerned. The bandgap shrinks remarkably as tensile strain increases, concomitantly with both electron and hole effective masses lowering, except for those of MoTe2. By contrast, when a compressive strain is applied, the bandgap and electron effective masses evolve at a much slower rate. In the meantime, the hole effective masses first increase and then decrease as strain goes from 0% to −4%. Nevertheless, the strain-induced bandgap shrinkage shall be strong enough to alter the strain-free bandgap energy sequence following MTe2, MSe2, MS2. Irrespective of the compound and strain, the absorption threshold and refraction peak occur at a lower energy in Figure 9. Results of h|BD| M-X +k|BD| X-X vs. ε 1 (0) of the MX 2 bilayer structures (a) under biaxial strain range of [−4%; +4%] and (b) under no strain at their maximal fitting coefficient of determination R 2 . Inset in (a) depicts the first-, secondand third-order fittings of R 2 vs. h/k. Inset in (b) depicts the third-order fitting of R 2 vs. h/k ratio. Arrows in between symbols "−" and "+" in (a) are used to represent the evolution from the biaxial compressive (−4%, −2%) strains to the tensile (+2%, +4%) ones.

Conclusions
The influence of biaxial strain on optoelectronic properties and electron density of transition metal dichalcogenides MX 2 (M = Mo, W and X = S, Se, Te) bilayers were thoroughly investigated for the first time using DFT calculations. When subjected to a strain going from a compressive (−4%, −2%) to tensile (+2%, +4%) one, the 2D materials' band structures and their corresponding effective masses of charge carriers undergo non-symmetrical changes as compression and tension are concerned. The bandgap shrinks remarkably as tensile strain increases, concomitantly with both electron and hole effective masses lowering, except for those of MoTe 2 . By contrast, when a compressive strain is applied, the bandgap and electron effective masses evolve at a much slower rate. In the meantime, the hole effective masses first increase and then decrease as strain goes from 0% to −4%. Nevertheless, the strain-induced bandgap shrinkage shall be strong enough to alter the strain-free bandgap energy sequence following MTe 2 , MSe 2 , MS 2 . Irrespective of the compound and strain, the absorption threshold and refraction peak occur at a lower energy in the in-plane (xx) direction than in the out-of-plane (zz) one. A strain effect is only visible in the xx-direction. For a given compound, the absorption threshold is lowered when subjected to a tensile strain. Irrespective of the metal atom and strain, the refraction peak intensity increases, and the absorption edge value decreases as the chalcogen atomic number increases. The strain effect on absorption and refraction is direction dependent. More precisely, these values are further improved in the xxand zz-direction when the compound undergoes tensile and compressive strains, respectively. From the determination of the bond degree (BD) at the bond critical points using QTAIM, it was found that the X-X BD decreases when the chalcogen atomic number increases and is nearly constant with applied strains. The kinetic energy per electron G/ρ at the bond critical point was also estimated. It is almost independent of strains with the lowest values for X = Te, which can be related to the highest polarizability evidenced in MoTe 2 from static relative permittivity values. A cubic relationship between the absolute BD summation of M-X and X-X bonds and the static relative permittivity was observed both in strained and unstrained MX 2 bilayers. After applying strain, the preponderant contribution of the X-X bonds in this relation under no strain was substantially reinforced, yielding the leading role to the X-X bonds instead of the M-X ones in the photovoltaic response. As the application of a biaxial strain on the chalcogenides bilayers allows for reproducing the effect of epitaxial strain on the absorber layer in a real device, these results can be valuable for the building of photovoltaic devices.