Band Nesting in Two-Dimensional Crystals: An Exceptionally Sensitive Probe of Strain

Band nesting occurs when conduction and valence bands are approximately equispaced over regions in the Brillouin zone. In two-dimensional materials, band nesting results in singularities of the joint density of states and thus in a strongly enhanced optical response at resonant frequencies. We exploit the high sensitivity of such resonances to small changes in the band structure to sensitively probe strain in semiconducting transition metal dichalcogenides (TMDs). We measure and calculate the polarization-resolved optical second harmonic generation (SHG) at the band nesting energies and present the first measurements of the energy-dependent nonlinear photoelastic effect in atomically thin TMDs (MoS2, MoSe2, WS2, and WSe2) combined with a theoretical analysis of the underlying processes. Experiment and theory are found to be in good qualitative agreement displaying a strong energy dependence of the SHG, which can be exploited to achieve exceptionally strong modulation of the SHG under strain. We attribute this sensitivity to a redistribution of the joint density of states for the optical response in the band nesting region. We predict that this exceptional strain sensitivity is a general property of all 2D materials with band nesting.

Data for WSe 2 and WS 2  Figure S1: Experimental (left) and theoretical (right) total SHG intensity I tot = π(A 2 + B 2 ) and distortion D = (A + B) 2 /(A − B) 2 for WS 2 (top) and WSe 2 (bottom) at measured and calculated strain values (see legend). Solid vertical lines mark the resonance frequency ω R = Ω R /2. D is plotted in the vicinity of ω R in scaled units (shaded area in theory shows corresponding energy range, for the experimental data the energy scale and peak width scale coincide). The insets show one typical angular resolved SHG signal at an energy where strain has a large effect (marked by dashed vertical line), with ϕ = 0 fixed to the direction of strain (see SI for data at additional energies). Different lattice orientations w.r.t. the direction of strain result in rotated and deformed angular resolved patterns.

S-2
Extended Data for MoSe 2 and MoS 2 We extend the data presented in the main text with additional experimental and theoretical data of the polarization resolved second harmonic generation at different energies.
The SHG signal in MoSe 2 under strain is distorted such that the lobes of the polarization resolved signal are smaller in the direction of strain ϕ = 0 than in the strain-orthogonal direction, corresponding to D < 1. The inverse behavior is observed in MoS 2 where the lobes in the direction of strain are larger than orthogonal to it, corresponding to D > 1.
Clearly, there are differences between theory and experiment. This is not surprising considering that we neglected many-body interactions beyond the PBE functional in DFT in the calculations. But we stress that the general trend is represented well already within the independent particle approximation, showing that band-nesting, which can be described at this level of theory, is already enough to explain the general trend in the data.
S-3  Divergence of optical response due to band nesting We give a brief review on band nesting and on the effect of dimensionality on the divergence of the joint density of states and the resulting features in the optical constants following reference. S1 The imaginary part of the dielectric function in a gapped material is given by where D is the dimension of the material under investigation, e is the elementary charge, m the electron mass,ê the direction of the laser polarization, ω the photon energy, C,V the energy of the conduction (C) and valence (V ) states and M C,V (k) the transition matrix element between them. If we assume M C,V (k) to be slowly varying over k we find that the strength of an optical transition is determined by the so-called joint density of states where d D S is a D-dimensional surface element on the surface defined by C (k) − V (k) =hω.
From Eq. S2 we see that special care has to be taken for points where vanishes. The surface element d D S counteracts divergences in the denominator. From this we can deduce that the effect of band nesting will be larger for lower dimensions. Using a Taylor expansion in k around the singularity in the denominator, Bassani and Parravicini S1 show that in band nesting regions the joint density of states only shows a kink in three dimensions and may diverge in two dimensions (at saddle points) and always diverges in one dimension.
This behavior has been shown to persist also for higher-order optical constants. S2,S3 S-5

Computational details
Linear response with and without excitonic corrections To elucidate the role of electron-hole interactions on the optical absorption, we compare the linear response spectrum in the independent particle approximation (IPA) to the one calculated by solving the Bethe-Salpeter-equation (BSE). In the BSE, we introduce a scissor operator which rigidly shifts the conduction band levels of the IPA by ∆ ≈ 1 eV to correct for the underestimation of the quasi-particle band gap in DFT. This upwards shift of the bands is counteracted by the electron-hole binding energy in the BSE. BSE Figure S3: Comparison of the linear response spectrum calculated in the independent particle approximation applied in this work and obtained by solving the Bethe-Salpeter equation.
As expected, the absorption spectrum is significantly altered close to the band gap as oscillator strength is shifted to two pronounced peaks, commonly referred to as A-and Bexciton [ Fig. S3]. Those two peaks stem from two energetically well separated eigenstates of the BSE and their energy is well below the quasi-particle band gap. In contrast to those welldefined low-energy excitons, excitonic corrections in the band nesting region have a different effect. As opposed to what the terminology C-exciton suggests, the band nesting peak results S-6 from a dense sequence of eigenstates of the BSE. S4 While those eigenstates show a significant energy shift w.r.t. the quasi-particle band gap, the shape of the resulting resonance in the IPA and the BSE look quite similar, further corroborating our approximation of neglecting explicit excitonic corrections in the second-order response.
Uniaxial strain and structural relaxation strain a (armchair) b (zig-zag) Poisson ration ν Figure S4: Schematics of the implementation of strain in the simulations. The cell is strained in a-direction and reduced in b-direction according to the Poisson ratio.
The simulations of the ground state atomic structure under strain are performed using the Vienna Ab-initio Simulation Package (VASP). S5,S6 Uniaxial strain is applied in the armchair direction (a-direction) of the crystal lattice (Fig. S4). In a first step we determine the Poisson ratio which is given by the relaxation of the unit cell orthogonal to the applied We then simulate strains between 0%-1% with the unit cell dimensions given by the determined Poisson ratio. With the unit cell in a/b-direction fixed, the atoms are then relaxed. A summary of the convergence parameters and the obtained lattice constants and Poisson ratios are given in Table S1.
S-7 Second order susceptibility χ The second order susceptibility is calculated according to Eq. (1) in the main text using the exciting code. S7 We include a scissor operator fitting the theoretically obtained band gap to the experimental one (Table 3 in S8 ). The convergence parameters used for the simulations are summarized in Table S2. Calculations on MoS 2 with the strain axis along the b-direction confirm that the same strain dependence is seen in both approaches.
ijk is a complex quantity where the real and the imaginary part are connected by a Kramers-Kronig relation. We have checked that for the energy range presented, all components of χ

Connecting experiment and theory
To connect the A-and B-parameters obtained experimentally and the calculated χ (2) ijk we compare the parallel polarization P || (ϕ) in both frameworks P || (ϕ) = |E| 2 · (A cos(3ϕ) + B cos(ϕ)) (3) where double subscripts are summed over. For the total distortion and intensity (Fig. S1), the parameters A and B are obtained from calculating P || (ϕ) according to Equ. S4 and then fitting Equ. S3.
For the k-resolved analysis we have to take a different approach. Using the mirror symmetry of the crystal which prevails in the simulational setup, we find that three distinct, non-zero matrix elements of χ xyy and χ (2) yyx . By comparing the analytical expressions in Equ. S4 and Equ. S3 we can extract a connection between A, B and χ Those relations were used to plot the distortion and the total intensity in k-space,

Analysis of the origin of the strain dependence
To investigate the origin of the strain dependence we performed restricted summations in k-space centered around different special points (Γ, K, K', Q) in the BZ and compared the qualitative behavior of the χ ijk -curves with strain. For MoS 2 , WS 2 and WSe 2 , we find that integrating around Γ already yields converged χ ijk -curves while for MoSe 2 no such region could be found. Instead, it was necessary to include large areas of the BZ including the Γ and Q points.

Linearization
We connect our findings to the recently proposed method for strain measurements in TMDs using SHG. S9 The concept relies on two assumptions: (i) SHG is sensitive to strain and (ii) there is a linear parametrization of A and B with strain. Our results highlight that sensitivity can be optimized by tuning the excitation wavelengths. However, considering the complex evolution of the band nesting region with strain does not seem to suggest a linear behavior. We assess the validity of a linearization by fitting a linear regression model to the (potentially non-linear) evolution of A and B with strain. The coefficient of determination S10   Figure S10: Coefficient of determination R 2 comparing a linear (purple circles) or a quadratic (green triangles) fit to the strain dependence of A (bottom panels) and B (top panels) as function of energy for all four investigated TMDs. The intensity is plotted alongside as a reference (dashed line).

Experimental SHG Setup
We used a tunable (1.3 eV to 1.6 eV) Ti:sapphire laser (200 fs, 76MHz) as fundamental for the SHG measurements. With a quarter wave plate the fundamental is circularly polarized which is optimized with great care using a polarimeter. This circularly polarized light is then linearly polarized under an angle φ with a broadband wire grid polarizer mounted on a rotation stage. The light is focused on the sample with a 100× objective lens and the second harmonic is collected in backscattering geometry with the same objective lens. Before the SHG is detected by an avalance photodiode (APD), it is passed through the same wire grid polarizer the fundamental has passed through. After the SHG light is detected, it is separated from the fundamental by a dichroic mirror and lowpass filters. is important for the homogeneity of the spatial strain distribution over the flake. We apply uniaxial strain to the TMD monolayer by three-point bending.

S-17
The dielectric environment of our sample is given by the relative dielectric constant of SU-8, which is ∼ 2.8. S11 Absorption measurements for such a dielectric environment for all four investigated materials can be found in the literature. S12,S13 Comparing TMD monolayers on bare silicon (with a dielectric constant of about 11.7 S14 ) to monolayers suspended in air shows that the dielectric environment only slightly affects the absorption features. S15