Directional dependence of band gap modulation via uniaxial strain in MoS2 and TiS3

Strain is widely employed to modulate the band structures of two-dimensional (2D) van der Waals (vdW) materials. Such band engineering with strain applied along different crystallographic directions, however, is less explored. Here, we investigate the band gap modulation of layered chalcogenides, MoS2 and TiS3, and the dependence of their band gaps on the directions of applied strain, using first-principles calculations. The band gap transition in MoS2 is found to reduce in energy linearly as a function of increasing tensile strain, with a weakly directional-dependent gradient, varying by 4.6 meV/% (from −52.7 ± 0.6 to −57.3 ± 0.1 meV/%) from the zigzag to armchair directions. Conversely, the band gap in TiS3 decreases with strain applied along the a lattice vector, but increases with strain applied in the perpendicular direction, with a non-linear strain-band gap relationship found between these limits. Analysis of the structure of the materials and character of the band edge states under strain helps explain the origins of the stark differences between MoS2 and TiS3. Our results provide new insights for strain engineering in 2D materials and the use of the direction of applied strain as another degree of freedom.


Introduction
2D semiconductors have emerged as a fertile ground for research in recent years spurred on by both fundamental interest and potential for applications in photonics, catalysis and nanoelectronics [1][2][3][4][5].Such materials have a layered structure, with the shared characteristic of high in-plane strength and weak out-of-plane van der Waals interactions between layers, enabling easy separation of layers to form mono-to few-layer samples.This can often result in interesting optoelectronic properties such as high carrier mobility and photoluminescence [6].Strain engineering has emerged as a method of varying the band structure in such materials to tune both optical and electronic properties for different applications due to the high strains such materials can withstand [7].For example, local strain such as wrinkles has been shown to cause a 'funneling' effect, moving carriers towards strained areas-a phenomenon potentially useful for single photon emission [7].More generally, the method has been extensively shown to modify the band gap in TMDs and other van der Waals materials, enabling the emission wavelength to be tuned via strain [8,9].
One such 2D van der Waals material, MoS 2 , has been well studied due to its band gap character changing with number of layers.The change in band gap character in MoS 2 with reduced number of layers is a result of a quantum confinement effect, as the sulfur states (near the Γ point) extend spatially across layers whereas the molybdenum states (near the K point) do not, due to the position of molybdenum within the unit cell.The result is that reducing layer thickness only changes the energies of states at points within the Brillouin zone (BZ) corresponding to sulfur states, dchanging the band structure more significantly at the Γ point than at the K point and leading to a change in band gap position in reciprocal space [10].This change in nature and position in k space of the band gap is also known to occur when strain is applied to the system [11].
TiS 3 , is a layered material with many similarities to MoS 2 , such as its direct band gap and semiconducting nature.It has also gained interest due to its smaller bandgap, i.e. longer emission wavelength in the infrared range [12].MoS 2 and TiS 3 can both be exfoliated as single-or N-layer flakes, with the MoS 2 band gap transitioning from direct to indirect as the number of layers are increased, accompanied by a reduction in photoluminescence (PL) yield [10], while the TiS 3 band gap character remains stable for mono and fewlayer systems.This is an important property since it would allow higher thickness tolerances, making TiS 3 particularly attractive for device and engineering applications.The band gap in both materials can be further manipulated with strain, enabling control over the position of the peak wavelength of the PL spectrum [7,13,14].However, the effect of the strain direction on the band structures is not well studied: a question which could impact viability of real-world applications if crystal alignment to strain direction must be considered.Heterostructures involving multiple vdW materials could also be affected, as strain can be introduced in the heterostructure growth or formation process.The impact of the strain could then vary depending upon its direction within the layers [15].The directionality effect could also impact the field of twistronics, where angle is a highly critical parameter due to the large variation in supercell size with small twist angle [16].In this paper, we investigate this open question for both TiS 3 and MoS 2 using density functional theory.
Differences in symmetry between these two materials are likely to cause different directional responses to strain.2H-MoS 2 has a hexagonal symmetry, with its atoms arranged in a trigonal prismatic structure and space group P63/mmc in its bulk form, which then loses inversion symmetry in monolayer (ML) form, with the space group changing to P m 6 2 [17].TiS 3 has a monoclinic structure with a P2 1 /m space group, making it more anisotropic [18], however this can be approximated to an orthorhombic structure for the monolayer system, as the angle of the cell perpendicular to layers can be neglected for the monolayer case.The bulk structures for MoS 2 and TiS 3 (figures 1(a)-(b) and (d)-(e) respectively) consist of layers stacked in the c-direction with weak van der Waals forces dominating between layers.For MoS 2 , two transitions are important: the bulk indirect transition and the ML direct transition.The bulk band gap transition in reciprocal space is located between the Γ point and the halfway point between the Γ and K points (figure 1(c)), which then changes to a direct transition at the K point when in the ML form [10,19].The band gap transition is consistently direct and at the Γ point in ML and few-layer TiS 3 systems (figure 1(f)) [12].The MoS 2 bulk and ML band gaps have been measured experimentally as 1.29 eV and 1.9 eV respectively using a combination of PL, optical absorption and photo-conductivity techniques [20], with density functional theory (DFT) predictions varying significantly depending on the approximations employed (discussed further below).
Reduction in the band gap energy with increased lattice strain in MoS 2 has been reported experimentally using different methods of introducing strain, many of which may also be applicable to TiS 3 [21][22][23], for example Manzeli et al applied strain using an atomic force microscope (AFM) tip and recorded a change in band gap measured via the piezoresistive effect [22].This measurement has been supplemented with first principles predictions at various levels of theory [24][25][26].Asymmetrical biaxial strain in MoS 2 involving some directional effects has been investigated computationally by Deng et al: fixing both lattice parameters a and b at various strain values (where strain along the two directions is not necessarily equal) to find the band gap in the a−b parameter space [27].
In this study, to model uniaxial strain we instead fix the unit cell for both materials along the direction of applied strain only, allowing for full structural optimisation in the perpendicular direction.Strain oriented along the two directions has been investigated for TiS 3 , however not for intermediate directions [14].Once the directional dependence of strain on the band gap is well understood, this determines if the alignment of a crystal relative to the strain axis necessary for the desired electrical and optical properties and performance in mechanically flexible devices, such as sensors, LEDs, transistors, and photodetectors [28,29].

Method
All DFT calculations were undertaken using a plane wave basis set, as implemented in the Vienna ab initio simulation package (VASP) [32].A gamma centred Monkhorst-Packgrid of 8 × 8 × 1 k points was used to sample the Brillouin zone for all calculations, with optimized structures converged to a force tolerance of 0.01 eV Å −1 .Both Perdew-Burke-Ernzerhof (PBE) and Heyd-Scuseria-Ernzerhof (HSE06) functionals were employed, with and without Grimme's D3 corrections (shown as +D3 after functional name) [31].For hybrid calculations, all supercells were first pre-optimised using PBE to allow for quicker convergence.Ideal values of D3 parameters for HSE06 are still an open area of research and could therefore not be sourced from literature, so this study used parameters quoted for the related hybrid functional PBE0 for calculations.Band calculations were then run on the optimized structures along a path through k space generated by SUMO package [33].Version 5.2 PBE plane wave potentials (PAW) were used for all calculations, with plane wave cutoff value of 520 eV.Spinorbit coupling was included for both geometry optimisation and band structure calculations.
Numerous DFT calculations have also been undertaken for MoS 2 using various levels of theory in order to calculate different properties.In plane lattice constants are generally close to the experimental values using either the generalised gradient approximation (GGA) or hybrid levels of theory, for example a = 3.17 Å (PBE) [36] or a = 3.18 Å (HSE06) [37,38].There is some variation in the literature relating to different code implementations and parameters.The out-ofplane lattice constant is governed by the weaker inter-layer coupling, so is generally overestimated if Grimme's dispersion corrections are not included [31,36,39].With regard to the choice of functional, one issue is that, while yielding reasonably accurate lattice structures, GGA functionals such as PBE are known to underestimate band gaps [40].This error can be reduced by using a hybrid functional which combines a GGA functional with a fraction of Hartree-Fock exchange [41].The lattice structure also affects the band gap, so for a layered material inclusion of an approximate treatment of van der Waals forces is necessary to enable more accurate lattice structure prediction.This is irrelevant in the monolayer systems, but is included for bulk systems [19,42].The exciton binding energy can lead to a large difference between the calculated electronic band gap and the optical band gap measured experimentally [43,44].GW calculations can predict a more accurate band structure when including the exciton binding energy, although trends in band gap should remain accurate for GGA and hybrid methods.The lattice parameters and band structure for TiS 3 were not calculated at different levels of theory to compare, as it was assumed that the materials are similar enough that the same functional for both materials would suffice.The validity of this assumption is further evidenced by the accuracy of lattice constants calculated to within 0.04 Å and 0.02 Å for a and b lattice constants respectively compared to experimental results [45].This also ensured the calculations could be more easily compared between the materials.
Optimisation of materials subject to uniaxial strain in the x direction was performed by a small modification to the VASP source code to ensure the stress matrix was only nonzero along the direction perpendicular to the applied strain, therefore only enabling optimization along this direction.The desired strain is then imposed by multiplying the components of the lattice vectors along the strain direction by the relevant strain, as shown in figure 2. Strain can be applied in arbitrary directions by rotating the supercell in the a-b plane such that the desired strain direction is aligned along the x direction.We further verified the method reaches an energy minimum by calculating the energy of structures with the angle between a and b vectors changed by ±0.5 degrees.
In order to gain further insight into the effect of strain on the band structure, band decomposed charge densities were calculated using the k points corresponding to the band edges and charge isosurfaces visualised in the VESTA package [30].Due to limitations with the k point grid sampling, it was not possible to use k points corresponding to the exact conduction band minimum so the closest point to the band edges in the k point grid was used for this type of analysis.

Results
The predictions of different functionals for lattice parameters and band gaps of bulk 2H-MoS 2 are compared in table 1.The in plane lattice constant, a, is close to that measured in experiment for all functionals tested.The out of plane lattice constant, c, was found to be far higher than experiment for both PBE and HSE06 unless D3 corrections were applied, which significantly reduced c to be close to the experimental value.This suggests PBE+D3 or HSE06+D3 would be good choices for this material.Table 1, The dependencies of direct and indirect band gaps on functional choice are also shown in table 1.Firstly, one can see there is a significant difference between the direct band gaps at the PBE and PBE with no spin-orbit coupling (SOC) levels of theory.Therefore, we include SOC in all subsequent calculations.It can also be seen from the difference between PBE and PBE+D3 functionals that the effect of D3 parameters was to reduce the indirect band gap, maintaining a similar direct band gap between the two cases.This was then mirrored in the hybrid calculations.Another effect of the inclusion of D3 parameters for both HSE06 and PBE functionals was a change in location of the conduction band minimum (CBM) from the K point to the K-Γ midpoint, bringing the band structure closer to profiles calculated in the literature [10].Finally, the effect of changing from PBE and HSE06 was to increase both the indirect and direct band gaps by a similar amount.Thus we conclude that the closest band gap to experiment is obtained using the HSE06+D3 functional.Although this approach predicts a fundamental gap slightly higher than the experimental band gap, the latter is based on optical measurement of the band gap, which is expected to be lower than the fundamental band gap by the exciton binding energy.
We then investigate the strain-engineered properties of ML MoS 2 and TiS 3 .The unit cells used for the ML were generated using the HSE06+D3 optimised bulk structures (figure 1) with a vacuum gap of 20 Å introduced, sufficiently large to remove interactions between periodic images.A symmetrically nonequivalent range of strain directions were sampled according to the unit cell symmetry.A 30 degree range and 90 degree range for MoS 2 and TiS 3 respectively were spanned (direction defined in degrees anticlockwise from the x direction shown in figure 2) and sampled every 15 degrees with various levels of uniaxial strain applied up to 5%.A larger range was necessary for TiS 3 due to the reduced symmetry in comparison to MoS 2 .Both direct and indirect transitions were investigated for MoS 2 to assess which transition has a lower energy, but this was unnecessary for TiS 3 as the direct transition has the lowest energy for both ML and few layer systems [12].
MoS 2 has two independent bonds per unit cell, one of which is increased in length while the other decreases with applied strain, as expected to maintain a positive poisson ratio.For each of these bonds, the difference in bond length Method of applying uniaxial strain by rotating the unit cell to align with the x direction then increasing the components of the lattice vectors along this direction (x and y directions shown in axes).Only y components of lattice vectors are allowed to relax during cell optimisation.MoS 2 and TiS 3 are both shown with strain directions defined as angle from the a lattice vector.Exaggerated strain is applied to (a) to illustrate approach.for 5% strain along 0degand 30degdirections is only 0.01 Å, suggesting a highly isotropic effect.TiS 3 has many more independent bonds per unit cell due to its more complex structure, and here we focus on the two most significant bonds for band gap determination (the bonds with highest band edge charge density localised along these two bonds, which are circled in blue and black in figure 4(b).These are oriented along two perpendicular directions, hence can be strained independently from each other.The change in length of these bonds for 5% strain along 0°and 90°is 0.2 Å and 0.1 Å for the Ti-Ti and Ti-S bonds respectively-an order of magnitude larger than the changes seen between strain directions for MoS 2 .Figure 3 shows how band gaps vary with strain applied in different directions in the ML systems.Band gaps in MoS 2 vary linearly, independent of direction.Least square linear fits were applied to direct and indirect gap trends separately and gradients with statistical error calculated from the fit.The strain at which the band gap changed from direct to indirect for each direction was found to be around 5.3%-5.4% for all three directions (although this is quite sensitive to the level of theory employed, so it is expected that the gradient of the band gap energy with respect to strain would be more useful for comparison to experiment [48]).In TiS 3 the story is different, as the band gap variation can only be reasonably approximated as linear with respect to strain for the extrema at 0 and 90 degrees, hence the gradient of the response is not a good figure of merit for this system.Interestingly, the band gap energy change (ΔE g ) with increasing strain changes from negative, when applied along the a lattice vector, to positive, when applied at 90 degrees to this, highlighting the anisotropic nature of the effect in TiS 3 in contrast to the MoS 2 case.The total energies per unit cell as a function of strain are also shown for both MoS 2 and TiS 3 in figure 3(c) and figure 3(d).It can be seen that there is little difference in energy between directions at the maximum strains for MoS 2 (13.2 meV).This is in contrast to the case of TiS 3 , where a far larger difference in energy 77 meV can be seen between 0 and 90 degrees.This suggests TiS 3 has a greater anisotropy in its elastic constants than is the case for MoS 2 .
Changes in structure in the two materials under strain can also be noted.In MoS 2 , there are two bonds per unit cell, both of similar length.Under five percent strain, one bond changes from 2.387 Å (unstrained) to 2.416 Å and 2.433 Å with strain along 0°and 30°directions respectively-corresponding to percentage changes in length of 1.2% and 2.0%.The other bond changes from 2.387 Å (unstrained) to 2.370 Å and 2.375 Å with strain along 0°and 30°directions respectivelycorresponding to percentage changes in length of −0.7% and −0.2%.
In TiS 3 , two relevant bonds are identified, as shown in figure 4. The 'Short bond' changes under 5 percent strain from 3.383 Å (unstrained) to 3.365 Å and 3.552 Å with strain along the 0°and 90°directions respectively-corresponding to percentage changes in length of 0.5% and 5%.The Long bond' changes from 2.627 Å (unstrained) to 2.725 Å and 2.617 Å with strain along the 0°and 90°directions respectively-with percentage changes in length of 3.7% and −0.4%.
The differences in structure with strain aligned along different directions can be seen to be more pronounced in TiS 3 compared to MoS 2 due to notable differences in symmetry.
To provide further insight into the weak dependence of the band gap on strain direction for MoS 2 , we analyse band decomposed charge densities corresponding to electronic states at the band edges (figure 4).At both the CBM and valence band maximum (VBM) in the unstrained system, the MoS 2 charge density is clearly not isotropic, resulting in some directional dependence on the band decomposed charge densities.However, the band decomposed charge density is localised mostly around the molybdenum atoms for both bands.This means that there is similar charge density along each of the bonds, suggesting straining each bond has a similar effect on the VBM/CBM energy.After straining, the 0°and 30°cases show similar charge densities, with only a small difference in the distribution in the CBM.This small difference explains the weak directional dependence of strain on the band gap.
The band decomposed charge density calculation results shown in figure 4 provide some explanation as to the different response to strain in TiS 3 compared with MoS 2 .The valence band shown in figure 4(b) has charge localised along a 'long' Ti-S bond oriented parallel to the a lattice vector (circled in black), while the conduction band has charge localised along a 'short' Ti-Ti bond perpendicular to this direction (circled in blue).This causes the charge density along these 'long' and 'short' bonds to be modified more when strain is oriented along the a and b lattice vectors respectively.If stretching both bonds causes an increase in band energy, the valence band energy will be increased with strain along the a lattice vector, while the conduction band will be increased with strain along the b lattice vector.This will result in a decrease in band gap for the former case and an increase for the latter.The origin for the non-linearity observed at intermediate directions is not immediately clear and is likely due to the complex interplay between these two effects when there is a strain component along both a and b directions.Only the high symmetry directions were included here to illustrate the directional variation, as these directions showed the largest differences in the character of band edge states.

Discussion
The initial hypothesis that the strain-band gap relationship would be close to isotropic for MoS 2 and highly anisotropic for TiS 3 was proven to be correct.The nonlinear transition between the high symmetry directions in TiS 3 was unexpected however.This builds upon work by Kang et al where only the extreme directions were investigated in TiS 3 , going further to compare to MoS 2 and provide insight as to how the differences between the two arise [49].A study by Deng et al also found a similar directional symmetry of band gap variation for biaxial strain in MoS 2 [50].
The inclusion of excitonic effects using the Bethe-Salpeter equation could enable PL spectra to be calculated for easier comparison to experimental results and the exciton binding energy to be included in calculations.However, these types of calculations were not used in this study due to their prohibitive computational cost [51].We expect the hybrid DFT approach used in this study remain reliable for studying the trends in the variation of band gaps with strain, although quantitative evaluation of energies is likely to have 100s of meVs offset (exciton binding energy) from experimental values.
Weak directional dependence of the strain on band structure could be useful in all application of uniaxial strain.Thankfully in the case of MoS 2 , we would not need to consider the crystal alignment to gain close to the maximum change in band structure.This is an advantage experimentally, as crystal alignment is likely to limit real-world applications, as it can be difficult and is often not easily reproducible.As other 2H phase transition metal dichalcogenides (TMDs) have similar hexagonal structures to MoS 2 , this weak directional is likely also extendable to these materials (including WS 2 , MoSe 2 , MoTe 2 and WSe 2 ).
The advantage of TiS 3 is the freedom to increase or decrease the band gap by rotation of the strain axis.It is not clear what level of maximum strain could be supported by TiS 3 experimentally along the b lattice vector, as the material often breaks into long 1D flakes.This could be investigated further experimentally to discover what level of band gap change can be induced in reality.This behaviour is likely extendable to similarly structured trichalcogenides such as ZrS 3 , ZrSe 3 and HfS 3 due to the structural origin of such behaviour.Experiments could be devised to investigate these predicted effects e.g. by applying strain along different directions in the crystal and monitoring their photoelectric properties.

Conclusion
The variation of band gap with uniaxial strain applied along different directions for both MoS 2 and TiS 3 was investigated by DFT calculation.MoS 2 was found to exhibit negligible anisotropy in band gap modulation, but strong anisotropy is predicted for TiS 3 , resulting in a change in the gradient of the band gap with respect to strain from negative to positive, as the strain direction is rotated from the a to the b lattice vector directions.A non linear response was observed for the intermediary directions in TiS 3 due to the complex interplay between the system asymmetry and strain, whereas the MoS 2 system maintained a linear behaviour for all strain directions, reflecting its more symmetric structure.From this we can conclude that the band gap in TiS 3 , and likely other trichalcogenides with similar structures, can be either increased or decreased by strain, whereas such effect is only limited to a small decrease in MoS 2 and other hexagonal dichalcogenides.Such directional dependent stain-induced bandgap engineering could be useful in optoelectronic applications, for example to achieve fine tuning of the wavelength in light emitting devices.Further work is necessary to identify the maximum level of strain that TiS 3 can be physically supported along these crystal directions.

Figure 1 .
Figure 1.Structures, band structures and brillouin zones of MoS 2 and TiS 3 (structures created using VESTA [30]).(a), (b) and (d), (e) show the structures for MoS 2 and TiS 3 respectively (Mo purple, S yellow, Ti blue).Two in-plane orientations are shown for TiS 3 to illustrate its anisotropic structure.Two layers are shown for each material to show the layered structure.(c) and (f) show the band structures with 2D brillouin zones shown in insets for MoS 2 and TiS 3 respectively calculated on 2H-MoS 2 at HSE06 level of theory with Grimme's D3 corrections including spin-orbit coupling [31].Direct (red arrows) and indirect (black arrows) transitions labelled.

Figure 2 .
Figure2.Method of applying uniaxial strain by rotating the unit cell to align with the x direction then increasing the components of the lattice vectors along this direction (x and y directions shown in axes).Only y components of lattice vectors are allowed to relax during cell optimisation.MoS 2 and TiS 3 are both shown with strain directions defined as angle from the a lattice vector.Exaggerated strain is applied to (a) to illustrate approach.

Figure 3 .
Figure 3. Changes in band gap transition energies compared to the unstrained system (ΔE g ) as a function of strain, with strain along different directions for (a) MoS 2 and (b) TiS 3 , and total energies for these strains and directions are shown for (c) MoS 2 and (d) TiS 3 .Modulations of the direct transition is shown for TiS 3 , and modulation of both direct and indirect transition shown for MoS 2 Linear fits are shown for (a) but lines are just a guide to the eye for (b)-(d).

Figure 4 .
Figure 4. Band decomposed charge densities in (a) MoS 2 and (b) TiS 3 for VBM and CBM, for the unstrained system, and 5% strain applied along high symmetry directions.Isosurfaces are visualised using VESTA [30] with isosurface levels of 2.5 × 10 −5 eÅ −3 for the VBM and 3 × 10 −5 eÅ −3 for the CBM.Relevant TiS 3 bonds (most affected by strain) for strain in different directions referred to as 'long bond' and 'short' bond are circled in black and blue respectively.The VBM and CBM charge densities for TiS 3 are shown viewed along the b and a lattice vectors respectively to highlight directions with largest changes (unnecessary for MoS 2 due to higher symmetry).

Table 1 .
Comparison of lattice parameters and band gap transitions predicted by different functionals for MoS 2 .E g Direct and E g Direct are the direct and indirect band gap transitions respectively.