# Interlayer exciton mediated second harmonic generation in bilayer MoS2

Shivangi Shree, Delphine Lagarde, Laurent Lombez, Cedric Robert, Andrea Balocchi, Kenji Watanabe, Takashi Taniguchi, Xavier Marie, Iann Gerber, Mikhail Glazov, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Shivangi Shree, Delphine Lagarde, Laurent Lombez, Cedric Robert, Andrea Balocchi, et al. Interlayer exciton mediated second harmonic generation in bilayer MoS2. Nature Communications, 2021, 12 (1), pp.6894. 10.1038/s41467-021-27213-8. hal-03863822

HAL Id: hal-03863822

https://hal.science/hal-03863822v1

Submitted on 27 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Interlayer exciton mediated second harmonic generation in bilayer MoS$_{2}$

Shivangi Shree¹, Delphine Lagarde¹, Laurent Lombez¹, Cedric Robert¹, Andrea

Balocchi¹, Kenji Watanabe², Takashi Taniguchi³, Xavier Marie¹, Iann C. Gerber¹,

Mikhail M. Glazov⁴,∗ Leonid E. Golub⁴, Bernhard Urbaszek¹,† and Ioannis Paradisanos¹

¹Université de Toulouse, INSA-CNRS-UPS, LPCNO, 135 Avenue Rangueil, 31077 Toulouse, France

²Research Center for Functional Materials, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan

³International Center for Materials Nanoarchitectonics,
National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan and

⁴Ioffe Institute, 194021 St. Petersburg, Russia

Second harmonic generation (SHG) is a non-linear optical process, where two photons coherently combine into one photon of twice their energy. Efficient SHG occurs for crystals with broken inversion symmetry, such as transition metal dichalcogenide monolayers. Here we show tuning of non-linear optical processes in an inversion symmetric crystal. This tunability is based on the unique properties of bilayer MoS₂, that shows strong optical oscillator strength for the intra- but also interlayer exciton resonances. As we tune the SHG signal onto these resonances by varying the laser energy, the SHG amplitude is enhanced by several orders of magnitude. In the resonant case the bilayer SHG signal reaches amplitudes comparable to the off-resonant signal from a monolayer. In applied electric fields the interlayer exciton energies can be tuned due to their in-built electric dipole via the Stark effect. As a result the interlayer exciton degeneracy is lifted and the bilayer SHG response is further enhanced by an additional two orders of magnitude, well reproduced by our model calculations.

**Introduction.**— Non-linear optical spectroscopy of atomically thin crystals gives crucial insights into their physical properties and investigates potential applications [1, 2]. Non-linear optics is based on the interaction of photons within the crystal. This can give rise, among other processes, to second harmonic generation (SHG) where two photons of identical energy combine into one photon of twice this energy [3, 4]. Currently existing applications in photonics and laser physics make use of this process in bulk crystals for example using potassium dihydrogen phosphate (KDP) or lithium niobate (LiNbO₃) [5].

SHG is vital both as a spectroscopic tool [6–8] and also for applications of a large class of materials such as II-VI and III-V semiconductors [9–11], nanotubes [12], magnetic- and non-magnetic layered materials [13–19], perovskites and antiferromagnetic oxides [20, 21].

The occurrence of SHG in a crystal is directly linked to its symmetry, where crystals with broken inversion symmetry can exhibit SHG and crystals which are inversion symmetric can typically not [5]. This dependence on symmetry makes non-linear optics in atomically thin crystals such as graphene a very rich field of research [22–28]. Although pristine graphene is an inversion symmetric crystal, the crystal symmetry or symmetry of the electronic states is potentially very sensitive to charges or single molecules in the layer vicinity [29], an imbalance in population of valley states [30], or electric currents [31]. However, graphene is a gapless material, so light-matter interaction cannot be controlled to the same extent as in atom-thin semiconductors based on transition metal dichalcogenides (TMD).

TMD monolayers such as MoS₂ and WSe₂ have broken

inversion symmetry and show intrinsically strong SHG signals [1, 32, 33]. Excitons are strongly bound in these materials [34] and govern optical processes also at room temperature, contrary to the model semiconductor GaAs [10]. This is crucial for non-linear optics as exciton resonances enhance light-matter interaction by several orders of magnitude [34]. As a result, in monolayers the excitonic contribution to SHG can be orders of magnitude higher than the intrinsic contribution from the crystal which exists also off-resonance [16, 17, 35–37]. An intrinsic challenge for exploiting SHG signals from nanostructures is the overall small signal, despite the giant response of a monolayer per unit thickness [33, 38, 39]. Hence, in the interest of both fundamental physics and applications we need to uncover the exact origins of SHG in layered semiconducting materials to enable further tuning and amplification.

In this work we demonstrate strong, tunable SHG based on exciton resonances in nominally inversion symmetric MoS₂ bilayers with 2H stacking: Varying the laser excitation energy allows us to address resonantly not only intralayer and but also interlayer exciton states with large oscillator strength [40], which makes MoS₂ bilayers particularly interesting. At these specific energies we generate SHG signals in MoS₂ bilayers which are (i) comparable in amplitude to off-resonance SHG monolayer signals and (ii) several orders of magnitude larger than reported in the literature for TMD bilayers with 2H stacking [32, 37, 41, 42]. The prominent role of interlayer excitons in MoS₂ bilayers allows for an additional SHG tuning mechanism: in applied electric fields we lift the degeneracy of the interlayer exciton states through the Stark effect. This results in an additional enhancement
<figure><img src="image_3.png" /><figcaption>FIG. 1: Second harmonic generation (SHG) spectroscopy in bilayer MoS$_2$ with 2H stacking in different dielectric environments. (a) Schematics of SHG in MoS$_2$ bilayers with 2H stacking. (b) SHG spectra plotted as a function of twice the laser energy $2E_L$. Each spectrum is a separate measurement and we plot all the SHG spectra for the different laser energies as in [35]. SHG experiments are performed on non-encapsulated 2H MoS$_2$ bilayers (laser power 4 mW) and hBN-encapsulated samples (10 mW), see panel (e). The main excitonic transitions are marked: intralayer A, B and interlayer IE. (c) Differential reflectivity spectrum of the same sample as a function of photon energy. (d) Polarization-resolved SHG intensity [32] for $2 \times E_L = 2.01$ eV in resonance with the interlayer exciton IE (black squares) energy of 2H MoS$_2$ bilayers. For comparison we show the laser reflection (red circles) which remains constant for each angle. (e,f) same as (b,c), but for bilayer samples encapsulated in hBN. Note that in panel (e) we have multiplied the SHG signal in the orange shaded region by 10, and in the yellow shaded region by 50 to discern all transitions on a linear scale.</figcaption></figure>

of the SHG signal at the interlayer exciton resonance by two orders of magnitude, surpassing the response of an ungated monolayer at this energy. We analyze the origin of the highly tunable SHG response through the linear polarization dependence of the SHG signal with respect to the crystallographic axis. We compare bilayer results with measurements in monolayers and trilayers located in the same van der Waals stacks. We present a model analysis of the emergence of the strong SHG in the bilayer, taking into account the impact of electric fields perpendicular to the bilayer as well as different sources of asymmetry.

**Results.**— Very recently details of the rich excitonic resonances in bilayer MoS$_2$ have been uncovered [40]. We aim to distinguish the possible SHG generation that comes from a weak, intrinsic, crystal related contribution from the possible excitonic contribution by using a tunable laser source. This allows comparing off-resonance with on-resonance experiments.

For the experimental results shown in Fig. 1 we investigate a 2H MoS$_2$ bilayer deposited on 90 nm of SiO$_2$ on top of a Si substrate, a very typical sample configuration. We use a pulsed Ti:Sapphire laser source coupled to an optical parametric oscillator (OPO) and we scan twice the laser energy (i.e. $2 \times E_L$) over an energy range
<figure><img src="image_4.png" /><figcaption>FIG. 2: Comparing SHG spectroscopy for MoS₂ mono-, bi- and trilayers. hBN encapsulated sample with staircase arrangement of MoS₂ mono-, bi- and trilayer. (a) SHG spectra plotted as a function of twice the laser excitation energy $2E_L$ for the monolayer. The main exciton transitions are marked, compare with Fig. 1e,f for the bilayer and the supplement for mono-and trilayers. (b) SHG spectra for bilayer, where intensity in orange shaded region is ×10, yellow shaded region ×50. (c) SHG spectra for trilayer. (d) SHG intensity compared for the same power density used on the same encapsulated sample for a laser energy of $E_L$ = 1.005 eV i.e. in resonance with IE in the bilayer, see panels (a-c), bright blue vertical line. (e) For a laser energy of $E_L$ = 0.95 eV (indicated by vertical bright grey line in panels a and b) we plot the emission spectra in the spectral vicinity of $2E_L$ for the monolayer (black) and for the bilayer (red) using the left and right y-axis, respectively.</figcaption></figure>

that corresponds to the main optical transitions in the MoS₂ bilayers, namely the interlayer and intralayer exciton resonances [43–45]. Experiments are performed at a temperature of $T = 4$ K in vacuum in a confocal microscope (excitation/detection spot diameter of the order of the wavelength), see Appendix. We plot in Fig. 1b a series of SHG spectra from the bilayer for a number of selected laser energies between $E_L$ = 0.95 eV to 1.13 eV in steps of ≈ 6 meV. Each SHG spectrum is a single peak that shifts in energy as we vary the excitation laser energy $E_L$. The spectral width of the SHG signal is limited by the laser pulse duration (ps).

We start with the off-resonance case: we do not detect any SHG signal within our typical integration time (2 minutes and excitation power of a few mW) when the laser energy ($E_L$) is tuned *below* half of the intralayer A:1s energy (i.e. the lowest lying direct optical transition

with large oscillator strength). But surprisingly for a crystal with inversion centre, as we increase the laser energy, we see clear resonances in the SHG response, i.e. a strong variation in SHG intensity as a function of laser energy $E_L$ in Fig. 1b. In Fig. 1c we plot the differential white light reflectivity spectrum of this sample on the same energy scale.

Comparing SHG spectra and white light reflectivity, we deduce that the maxima in the SHG experiments occur when $2 \times E_L$ is in resonance with the intralayer A and B exciton energies, separated by about 180 meV in 2H stacked MoS₂ bilayers [46]. For intralayer excitons the electron and the hole reside within the same layer. In the energy range between A and B intralayer excitons we also observe clearly that the SHG response is strongly enhanced by the interlayer exciton IE transitions [43–45]. In the absence of electric fields, there are two degenerate
<figure><img src="image_5.png" /><figcaption>FIG. 3: Tuning the SHG signal of a 2H MoS₂ bilayer in applied electric fields - Stark splitting. Sample with contacts to apply vertical electric field $F_z$ (a) SHG spectra from a gated bilayer, covering intralayer A and interlayer IE exciton resonances at $F_z = 0$ MV/cm (blue spectra) and $F_z = 0.17$ MV/cm (red spectra) as a function of twice the laser energy $2E_L$. (b) comparison of SHG signal for interlayer exciton IE (solid circles) versus intralayer exciton A (open triangle) as a function of applied electric field $F_z$. The SHG amplitude of IE is fitted with a quadratic function $\propto F_z^2$ (red line). (c) white light reflectivity for three different electric field values $F_z = 0, 0.07, 0.13$ MV/cm. Here we plot here the first derivative of the reflectivity spectra for a better estimation of the IE energies. (d) Schematic of a bilayer with 2 distinct interlayer configurations, that have opposite permanent electric dipoles $\pm\mu$. (e) SHG spectroscopy for $F_z = 0, 0.07, 0.13$ MV/cm tuning twice the laser energy $2 \times E_L$ across the interlayer exciton resonance IE.</figcaption></figure>

interlayer exciton configurations, with a hole delocalized over both layers and an electron in the top or bottom layer, respectively, as shown in Fig. 3d. Interlayer excitons in homobilayer MoS₂ have large oscillator strength and their energy can be tuned over a wide range by the Stark effect [40, 47]. We explore this tunability in experiments shown below (Fig. 3).

The surprisingly strong SHG signal, which occurs in our experiments only at the excitonic resonances of the bilayer, could suggest a symmetry breaking in our non-encapsulated sample. A possible reason can be a drastically different dielectric environment experienced by top and bottom MoS₂ layers: For the results in Fig. 1b the bottom MoS₂ layer is on SiO₂ and the top layer is nominally uncovered (we perform our experiments in vacuum). To verify the impact of the dielectric environment, we have fabricated MoS₂ bilayer samples encapsulated in high quality hexagonal BN [48]. This leads to a more symmetric dielectric environment [49–51], therefore eliminating potential sources of symmetry breaking. In

addition, we also exclude the possibility of adsorbates or molecules at the sample surface in encapsulated samples which might impact the measurements as in surface SHG studies [52] and suggested in graphene [23, 29]. In Fig. 1e we plot the SHG spectra as a function of twice the laser energy $2 \times E_L$ using a smaller step size of $\approx 3$ meV between two adjacent SHG spectra. The SHG signal coming from hBN is several orders of magnitude weaker than from the MoS₂ bilayer, see [32] and our data in the supplement. Also for this sample we do not detect SHG signals below 1.90 eV (intralayer A:1s resonance). As the laser energy is increased, we observe, as for the non-encapsulated sample, clear resonances. Coming back to our initial motivation, our experiments in Fig. 1e show that also hBN encapsulated MoS₂ bilayers in a far more symmetric dielectric environment show strong SHG at the exciton resonances. These resonances are spectrally narrower in the encapsulated samples, see Fig. 1e, as compared to the non-encapsulated samples (Fig. 1b). This is consistent with TMD exciton transitions in high
quality hBN encapsulated layers [53]. Comparison with white light reflectivity data in Fig. 1f allows attributing the different interlayer and intralayer transitions. The slight overall shift in energy of the exciton transitions for encapsulated versus non-encapsulated samples comes mainly from renormalization of all the Coulomb energies as the effective dielectric constant is different for the two bilayer samples. In addition to previously identified exciton transitions a new transition at 2.05 eV emerges in SHG spectroscopy, that we tentatively attribute to the A:2s state of the intralayer excitons. High field magnetoptics is needed to confirm this assignment [54] through measuring the diamagnetic shift.

In Fig. 1d we plot the polarization dependence of the SHG response at the interlayer exciton resonance. We excite with linearly polarized light and the strength of the SHG signal collected in the same polarization depends on how the crystallographic axes are aligned with respect to the laser polarization. We clearly observe a 6-fold rotational symmetry expected for the space group of a 2H bilayer [32], see discussion below. This polarization dependence is a strong indication that the SHG signal is due to intrinsic effects linked to crystal and exciton symmetry. We have performed measurements at additional exciton resonances that give the same polarization dependence. Although in principle defects with 6-fold symmetry can exist in the crystal, their direct role at the SHG at exciton resonances is unlikely to be dominant, as discussed below. We conclude that the measured SHG signal is not due to extrinsic effects such as impurities deposited at the surface, which can generate SHG signals unrelated to the crystal symmetry itself as studied in surface SHG [52, 55, 56], for graphene [23, 29] and hBN [15].

We provide in Fig. 2a-c comparison for hBN encapsulated mono-, bi- and trilayers, for the same sample used in [40]. Our aim is to compare the surprising SHG amplitude for bilayers with the monolayer and trilayers, for which crystal inversion symmetry is broken and as a consequence strong SHG is expected [32]. Strikingly, for all the experiments on the mono-, bi- and trilayers, we see that the SHG signal is orders of magnitudes enhanced when twice the laser energy $2 \times E_L$ is in resonance with an excitonic transition, as compared with a non-resonant situation.

Previous reports on SHG in MoS$_2$ bilayers did not focus on exciton resonances [1, 32, 37, 41, 42, 57] and hence signals from monolayers were 3 orders of magnitude higher than for bilayers in the off-resonant case. From our measurements in panels Fig. 2a and b we deduce that $E_L=0.95$ eV (i.e. SHG energy of 1.9 eV) is below the A-intralayer exciton resonance for monolayers and bilayers. We directly compare the measured SHG signal at this laser energy for monolayers and bilayers, see peaks in Fig. 2e, and find that the monolayer signal is indeed 2 orders of magnitude larger than the bilayer

signal. The situation is drastically different as we change laser energy: At $E_L=1.005$ eV twice the laser energy $2 \times E_L$ is in resonance with the interlayer IE transition of the bilayer resulting in a strong amplification of the SHG signal while being non-resonant for a monolayer, see comparison in Fig. 2d. These different situations with respect to the exciton resonance result in a comparable overall SHG amplitude for mono- and bilayers. Note that Ref. [41] discusses a broad resonance with the high energy C-exciton, but no analysis of intra- versus interlayer excitons is provided, as details of these highly tunable transitions were only revealed very recently.

The trilayer crystal has no inversion centre, as the monolayer, and SHG can be expected for all energies. In addition, there exist several intralayer and interlayer exciton resonances [40] that can potentially enhance SHG. For the measurements in Fig. 2c we find an SHG signal with strong amplitude variations over the investigated energy range 1.87 and 2.15 eV. We identify the two intralayer transitions for excitons in the central layer ($A_{L2}$) and in the outer layers ($A_{L1}$ and $A_{L3}$) respectively, between 1.9 and 1.92 eV. Local maxima around 1.98 eV can be attributed to interlayer excitons, with B-exciton contributions around 2.1 eV. The exciton transitions for mono- and trilayers are identified by comparing with white light reflectivity [40] shown in the supplement.

By choosing a suitable laser energy, we have shown that the bilayer SHG signal (on-resonance) can reach the same order of magnitude as in the monolayer (off-resonance). We now show that by applying an electric field normal to the bilayer, we can further increase the SHG signal in particular at the interlayer exciton resonance. The interlayer exciton has an in-built, static electric dipole and high oscillator strength, which makes it visible in absorption in Figs. 1c,f and 3c. This strong absorption feature is highly tunable in energy through the Stark effect [40]. In Fig. 3 we show that the SHG signal is tunable in amplitude and also spectrally when an external electric field is applied to the bilayer. In previous works the effect of the electric field application on SHG response of 2H bilayers [41, 42] has been addressed only in terms of the impact of doping and crystal symmetry breaking, while interlayer exciton tuning as in our resonant SHG experiment has not been observed.

In Fig. 3a we compare the SHG signal of a bilayer with and without an applied electric field. For this gated sample, the IE SHG signal reaches 5 counts/s at $F_z = 0$. As we apply a field of $F_z = 0.17$ MV/cm this signal increases by a factor of 25 to about 125 counts/s at the IE resonance maximum. To study this tuning in more detail, we plot in Fig. 3b the SHG amplitude as a function of the applied electric field for six separate experiments. We see a quadratic increase of the IE SHG signal as a function of applied electric field $F_z$. This quadratic increase is due to mixing with intralayer excitons, as we show below.

For the intralayer exciton we record a much weaker
increase of the SHG signal as $F_z$ is increased. At $F_z = 0.17$ MV/cm the SHG resonance shows a slight red-shift (Fig. 3a), possibly indicating charging effect, i.e., a shift towards trion transitions. Charging effects as a possible source of symmetry breaking in bilayer SHG response are discussed for WSe$_2$ [42], our data demonstrate that this has negligible effect on the SHG response at the inter-layer exciton energy.

So far we discussed the amplitude of the SHG response of the IE resonance as a function of the applied electric field $F_z$, we now turn our attention to the spectral range over which this amplification occurs. In a 2H homobilayer at $F_z = 0$ two degenerate IE configurations exist: with a hole delocalized over both layers bound to an electron either in the bottom (IE$_1$) or top (IE$_2$) layer, see sketch in Fig. 3d. Application of a non-zero $F_z$ lifts the degeneracy of IE$_1$ and IE$_2$ as their static, permanent dipoles point in opposite directions. This leads to a Stark shift to lower and higher energy, respectively [40]. This lifting of the degeneracy can be seen in Fig. 3c in white light reflectivity experiments. We extract a splitting between IE$_1$ and IE$_2$ of roughly 10 meV for $F_z = 0.13$ MV/cm. In Fig. 3e we plot SHG spectroscopy results at the IE resonance for $F_z = 0, 0.07$ and $0.13$ MV/cm, respectively. In addition to the increase in SHG amplitude with $F_z$, we also observe that the spectral range for which amplification is observed is broadened. This broadening is the consequence of the energy splitting between IE$_1$ and IE$_2$ transitions in an applied electric $F_z$. The broadening can also be seen directly in Fig. 3a. Tuning both amplitude and spectral range of the SHG response of the bilayer is a direct consequence of exploiting the properties of inter-layer excitons and will be further discussed below.

**Discussion.**— We first focus on the model considerations for SHG mediated by the interlayer excitons in bilayer MoS$_2$. Here the SHG signal is very sensitive to the electric field, being one of the most striking experimental observations. Further we briefly address other contributions to the SHG. The detailed theory is presented in the supplement.

A TMD monolayer is described by the noncentrosymmetric point group $D_{3h}$ which lacks an inversion center and allows for the SHG. An ideal homobilayer in absence of external fields has an inversion symmetry and is described by the $D_{3d}$ point group where the SHG is forbidden. Qualitatively, this is because the constituent monolayers in the bilayer are rotated by an angle of $\pi$ with respect to each other in 2H stacking [58], and the electric field at a double frequency generated in one monolayer is compensated by the contribution of the other, Fig. 1a.

An electric field $F$ perpendicular to the monolayers breaks the symmetry and enables the SHG described by a phenomenological relation

$$P_{2\omega,x} = \chi(E_{\omega,x}^2 - E_{\omega,y}^2), \quad P_{2\omega,y} = -2\chi E_{\omega,x} E_{\omega,y}, \quad (1)$$

where $E_\omega$ is the incident field and $P_{2\omega}$ is the polariza-

tion induced at a double frequency, $x$ and $y$ are the main in-plane axes of the structure. Since the bilayer is centrosymmetric, the non-linear susceptibility $\chi$ in Eq. (1) changes sign at the space inversion, particularly, $\chi(-F_z) = -\chi(F_z)$. Polarization dependence of the SHG predicted by Eq. (1) is observed experimentally, Fig. 1(d), showing 6-fold symmetry of the crystal structure.

Intralayer excitons show negligible Stark shifts in applied fields [40], as their static dipole moment is negligible, as for excitons in separate monolayers [59, 60]. Therefore, due to their strong in-built electric dipole the main effect of the field $F_z$ is on the indirect excitons $IE_1$ and $IE_2$ whose degeneracy is lifted yielding $E_{IE_i} = E_{IE} \mp F_z \mu$, where $\pm \mu$ is the static dipole moment of the corresponding IE, Fig. 3(d). As a result of the energy splitting (Fig. 3c), the SHG contributions from $IE_1$ and $IE_2$ do not cancel, resulting in a strong SHG signal as the energy separation between $IE_1$ and $IE_2$ increases (Fig. 3a). The second-order susceptibility in the vicinity of the IE resonances in the linear-in-$F_z$ regime can be calculated taking into account two-photon excitation and coherent one-photon emission from the exciton states [36] as (see supplement for details):

$$\chi_{\text{IE}} = -2\mu F_z \frac{|T|^2 C_2 |\Phi_{\text{B:1s}}(0)|^2}{\Delta^2 \left(2\hbar\omega - E_{\text{IE}} + \frac{|T|^2}{\Delta} + i\Gamma_{\text{IE}}\right)^2}. \quad (2)$$

Here $C_2$ is the combination of the interband momentum matrix element, electron charge, free electron mass and the band gap, $\Phi_{\text{B:1s}}(0)$ is the envelope function of the B:1s exciton at the coinciding electron and hole coordinates, $T$ is the hole tunneling matrix element between the layers, $\Delta = E_B - E_{\text{IE}}$ is the energy splitting between the IE and B excitons (which are mixed due to the hole tunneling [43, 45]), and $\Gamma_{\text{IE}}$ is the damping of the IE resonance. As demonstrated in previous works [40, 47], the IE exciton optical activity is mainly due to the mixing with B-excitons, thus, in Eq. (2) the ratio $|T/\Delta|^2$ accounts for the IE-B exciton mixing in the two-photon excitation and single-photon emission channel. Equation (2) holds provided $|\mu F_z| \ll \Gamma_{\text{IE}}$, otherwise two peaks in the SHG spectra split by $2|\mu F_z|$ are expected to split (see supplement for details).

Two important conclusions can be drawn from Eqs. (1) and (2): First, $P_{2\omega}$ is maximized when twice the laser energy is resonant with an indirect exciton state as $2 \times E_L = 2\hbar\omega = E_{\text{IE}}$. Second, the intensity of the SHG scales as $F_z^2$. Experimental SHG intensities in Fig. 3a,e and also in Figs. 1b,e,2b are strongly enhanced at the IE resonance. In the measurements in Fig. 3b the SHG amplitude grows quadratically as a function of the applied $F_z$. So both experimental observations strongly support our analysis.

Importantly, the experiment demonstrates also the SHG on other excitonic species in bilayers, including weakly field dependent effect on intralayer A-excitons,
Fig. 3(a,b), and also the SHG at zero field on intralayer A:2s and B-excitons. A weaker effect of the electric field *F* on the A-excitons can be readily understood taking into account the fact that the A-excitons are mixed, e.g., with the interlayer B-excitons via hole tunneling. Corresponding energy distance $\Delta' = E_{\text{IE(B)}} - E_A \gg \Delta$ resulting in a weaker susceptibility to $F_z$.

Interestingly, the SHG is also observed at the exciton resonances at $F_z = 0$, i.e., in the nominally centrosymmetric situation. Although the signals are smaller than those at the resonances in non-centrosymmetric monolayer and trilayer samples, they are comparable to the non-resonant background in monolayers and trilayers, cf. Fig. 2. Possible origin of SHG in this case could be the quadrupolar or magneto-dipolar SHG, where akin to graphene case [22, 26], the effect is related to the light wavevector $q_z$ at the normal incidence. The symmetry analysis again yields Eq. (1) with $\chi \propto q_z$. Particularly, a (tiny) phase difference of $\phi \sim q_z a$ where $a$ is the interlayer distance in the bilayer yields an imbalance of the contributions of individual monolayers and results in the resonant contributions at A and B excitons similar to calculated in [36] for monolayers, but being by a factor of $\phi \sim 10^{-2}...10^{-3}$ smaller (see supplement). Our estimates show, however, that this effect is too small to account for the surprising experimental findings which demonstrate that the resonant susceptibilities of the bilayer are roughly one order of magnitude smaller than those of the monolayer.

Thus, we arrive at the conclusion that, despite relatively symmetric environment of our bilayer sample, the structure lacks an inversion center. Possible options are (i) small built-in electric fields and (ii) inequivalence of the intralayer excitons in monolayers forming a bilayer. Option (i) results in the replacement of $F_z$ by $F_z + F_0$ where $F_0$ is the normal component of the build-in field. In this case at $F_z = -F_0$ the effectively symmetric situation can be realized. While the measured dependence of the SHG intensity for IE, Fig. 3, does not contradict this scenario, note the minimum at $F_z \approx -0.035$ MV/cm, the A:1s intensity does not significantly drop in this region. Alternatively, we may suppose that the build-in field is inhomogeneous in the sample plane within our detection spot and the contributions of the A:1s and IE states come from slightly different nanoscopic regions.

Option (ii) implies that the energies of A- and B-excitons or their oscillator strengths are inequivalent in the two constituent monolayers that form the bilayer. This may be due to the dielectric disorder, which, although suppressed, can still be present in state-of-the-art, hBN encapsulated samples [50]. In this situation the IE is mainly activated by the electric field, while intralayer excitons are less sensitive to $F_z$ (see Supplement). Inhomogeneous broadening effects due to disorder, impurities, etc. are inevitable [50], in particular in MoS$_2$ samples with different type of defects

reaching densities of $10^{10}$ to $10^{12}$ per cm$^{-2}$ [51]. But the crucial parameter is the ratio between the size of the inhomogeneity, $l$ and the laser spot diameter, $\delta$. For $l \ll \delta$, inhomogeneous broadening effects will effectively cancel out and the structure will effectively maintain the space inversion. However, for $l \gtrsim \delta$, the inhomogeneities will lift the inversion symmetry, thus a contribution of one monolayer will dominate over the other and SHG occurs. Sizeable effects due to disorder require a comparable or larger spatial scale than the excitation spot. Inhomogeneities on the order of one micrometer can occur even in hBN encapsulated samples [50]. It is important to note that sample inhomogeneity does not always translate in a large shift of the excitonic transition energy, as bandgap shifts and exciton binding energy changes partially compensate [50]. We stress that the six-fold symmetry of the SHG signal (see Supplement) clearly rules out, e.g., an in-plane asymmetry due to the strain or in-plane fields [19].

**Conclusion.**— In summary, we show strong and electrically tunable exciton-mediated SHG in 2H MoS$_2$ bilayers that can surpass the off-resonance monolayer signal. Drastic enhancement of the SHG amplitude is observed when twice the laser energy is in resonance with the excitonic transitions. At the interlayer exciton resonance, we tune the SHG signal by over an order of magnitude in electric fields applied perpendicular to the layer and demonstrate that the spectral width of the SHG resonance increases. We relate the SHG in applied electric field to the Stark splitting of the interlayer exciton and its mixing with intralayer excitons. Interlayer excitons form due to coupling between two layers. This coupling depends on several parameters such as layer distance and twist angle. Changing these parameters therefore allows tuning of SHG in van der Waals stacks. We identify several possible mechanisms for the strong zero-field SHG signal, such as small in-built electric fields in the structure. Disorder, possibly related to defects in the layer or in the dielectric environment, could result in non-equivalence of the intralayer excitons in the different layers within the optical detection spot on our sample structure. Our scheme for tuning SHG can be applied to a variety of other systems, such as homobilayer MoSe$_2$ with strong interlayer exciton resonances [61], hetero bilayers which host interlayer excitons with high oscillator strength for example MoSe$_2$/WSe$_2$ with hybridized conduction states [62] and WSe$_2$/WS$_2$ with hybridized valence states [63]. The impact of concentration and type of defects states on the non-linear optical properties merits further investigation [51, 64].
# METHODS

**Sample fabrication.-** 90 nm thick SiO₂/Si substrates were cleaned for 10 minutes in aceton and isopropanol using an ultrasonication bath and were subsequently exposed in oxygen plasma for 60 seconds. Bulk MoS₂ (2D Semiconductors) was first exfoliated on Nitto Denko tape and the exfoliated areas were attached on a polydimethylsiloxane (PDMS) stamp, supported by a microscope glass slide. Monolayers, 2H-bilayers and trilayers were identified based on the optical contrast under an optical microscope prior to transfer on the SiO₂/Si substrate. For the hBN encapsulated samples, hBN flakes were first exfoliated on a Nitto Denko tape from high quality bulk crystals [48] while the same PDMS-assisted transfer process on SiO₂/Si substrates was followed. A staircase sample of monolayers, 2H-bilayers and trilayers was subsequently transferred and capped in hBN. Between every transfer step, annealing at 150°C was applied for 60 min. For the electric field device, the same process was followed including the additional transfer of few-layered graphite (FLG) flakes. FLG flakes were aligned according to Au contacts on a pre-patterned substrate to act as electrodes for the electric field measurements. The precise sequence of the complete stack from bottom to top includes hBN/FLG/hBN/2H-MoS₂/hBN/FLG.

**Optical spectroscopy Set-up.-** Optical spectroscopy is performed in a home-built micro-spectroscopy set-up assembled around a closed-cycle, low vibration helium cryostat with a temperature controller (T = 4 K to 300 K). For SHG measurements we use ps pulses, generated by a tunable optical parametric oscillator (OPO) synchronously pumped by a mode-locked Ti:sapphire laser as in [17]. SHG signal is collected in reflection geometry. For low temperature white light reflectance measurements a white light source; a halogen lamp is used with a stabilized power, focused initially on a pin-hole that is imaged on the sample. The emitted and/or reflected light was dispersed in a spectrometer and detected by a Si-CCD camera. The excitation/detection spot diameter is 1 µm, i.e., smaller than the typical size of the homobilayers. We obtained differential reflectivity from reflectivity spectra as $(R_{ML} - R_{sub})/R_{sub}$, where $R_{ML}$ is the intensity reflection coefficient of the sample with the MoS₂ layer and $R_{sub}$ is the reflection coefficient of the hBN/SiO₂ stack.

# SUPPLEMENT

## Sample Details

Van der Waals heterostructures (vdWHs) were fabricated by stacking two-dimensional materials via a dry-

transfer technique. 2H-MoS₂ crystals (2D Semiconductors) and synthetic hBN [48] were subjected to micromechanical cleavage on Nitto Denko tape, then exfoliated again on a polydimethylsiloxane (PDMS) stamp placed on a glass slide for optical inspection. Sequential deterministic stamping of the selected flakes was then applied to build the complete stack. Optical images of hBN encapsulated but also bare MoS₂ in SiO₂ are shown in Fig. 4. The thickness of the bottom hBN layers was selected to optimize the oscillator strength of the interlayer exciton (IE) in bilayer MoS₂ [65].

## Experimental setup

The SHG spectra were recorded with the setup sketched in Fig. 5. An Optical Parametric Oscillator (OPO) is aligned to a closed-cycle cryostat to excite the sample under investigation. A combination of linear polarizers and halfwave plates allows the control of excitation and detection polarization for the polarization-resolved measurements. The light is focused onto the sample at 4.2 K using a microscope objective (NA=0.8). The position of the sample with respect to the focus can be adjusted with cryogenic nanopositionners. The reflected light from the sample is sent to a spectrometer with a 150 grooves per millimeter grating. The spectra are recorded by a liquid-nitrogen cooled charged coupled device (CCD) array. Low temperature reflectance measurements were performed using a halogen lamp as a white-light source with a stabilized power supply focused initially on a pin-hole that is imaged on the sample. The reflected light was dispersed in a spectrometer and detected by the same Si-CCD camera. The excitation/detection spot diameter is ≈ 1µm, i.e. smaller than the typical diameter of the sample.

## Assignment of resonances in monolayers and trilayers

Monolayers and trilayers have also been investigated via SHG spectroscopy. For the assignment of the different excitonic resonances in the SHG spectra, differential reflectance has also been collected from the same sample areas (compare top and bottom spectra in Fig. 6). Besides 1s states, excited states of the A exciton can be observed in SHG and verified by reflectivity. Note that interlayer exciton states are clear in the SHG spectra of trilayers, see also [40].

## Background SHG

For the gated bilayer device it is important to identify the background contribution in the SHG. Few-layered
<figure><img src="image_6.png" /><figcaption>FIG. 4: Optical images of hBN encapsulated and bare MoS₂ on SiO₂.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5: Experimental setup used for the collection of the SHG spectra.</figcaption></figure>

graphene (FLG), as well as hBN can provide a contribu-
tion in the SHG spectra. In Fig. 7 we show the noise level
of the CCD camera as well as the background SHG orig-
inating from hBN and FLG under zero and finite electric
fields. In this measurement, we collected the SHG com-
ing from the same stack of the gated sample, but without
MoS₂. The background SHG is weak (< 2counts/second)

compared to the signal from MoS₂ and it is wavelength
independent. This unambiguously confirms the spectral
signature of the excitonic resonances in SHG of MoS₂
bilayers.

Polarization-resolved SHG

Polarization-resolved SHG has been applied to verify
the symmetry and crystallographic orientation of MoS2
bilayers [32, 33]. The intensity of the SHG as a func-
tion of the polarization angle of the fundamental laser is
shown in Fig. 8. Notice that similar angle dependence
is acquired when the excitation is in resonance with the
A:1s and IE states. This SHG intensity is dependent on
the crystallographic orientation of the sample. The sym-
metric 6-fold co- (XX) and cross-linear (XY) polarization
over 2π reveals that any contribution from strain effects
and low symmetry defects is negligible [15, 19].

hBN nonlinear contribution

In the encapsulated bilayers with hBN, we measure
the SHG strength of the surrounding hBN in the total
signal. In Fig. 9 we compare the SHG intensity coming
from the surrounding hBN with the MoS₂ bilayer when
the excitation is kept constant and tuned into resonance
with the A:1s state. The SHG coming from MoS₂ is two
orders of magnitude stronger.

Validation of the nonlinear processes

Power dependent measurements have been performed
in MoS₂ bilayers. In Fig. 10 the SHG intensity is plotted
as a function of laser power. The excitation energy is
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6: Comparison between SHG and differential reflectivity spectra of monolayer and trilayer MoS₂.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 7: Background SHG contribution from the device excluding MoS$_2$ at zero and finite electric field.</figcaption></figure>

tuned in resonance with the A:1s state. The quadratic power law confirms the two-photon process.

## THEORY

This section contains the phenomenological and microscopic theory describing the SHG in bilayer structures.

### SHG in an ideal centrosymmetric bilayer

Bilayer has a $D_{3d}$ point symmetry with an inversion center. For the light propagating along the normal $z$ there is a contribution caused by the radiation vector $q_z$ (quadrupolar contribution) which results from the inhomogeneity of the incident and emitted electric field distribution:

$$P_{2\omega,x} = q_z\chi'(E_{\omega,x}^2 - E_{\omega,y}^2), \quad P_{2\omega,y} = -2q_z\chi'E_{\omega,x}E_{\omega,y}. \quad (3)$$

To illustrate the effect let us study the electrodynamic mechanism leading to Eq. (3). For this mechanism we consider both monolayers forming a bilayer as electronically decoupled, the distance between the layers is $a$. Disregarding multiple reflections of the second harmonics radiation from both monolayers we obtain the reflected and transmitted amplitudes of the second harmonics in the form

$$E_{2\omega}^r = S\chi E_{\omega}^2(1 - e^{i(q_{2\omega} + 2q_{\omega})a}), \quad (4a)$$

$$E_{2\omega}^t = S\chi E_{\omega}^2(e^{iq_{2\omega}a} - e^{i2q_{\omega}a}). \quad (4b)$$

Here $S$ is the coefficient describing the propagation of the SHG field from the sample, $\chi$ is the second-order susceptibility of a monolayer [36], and $q_\Omega = \Omega n_\Omega / c$ with $n_\Omega$ being the refraction index at the frequency $\Omega$. It is noteworthy that for a bilayer the products $q_\omega a, q_{2\omega} a \ll 1$. As a result, the intensity of the SHG from the bilayer is suppressed as compared with the monolayer by a factor of $\sim q_{2\omega} a \ll 1$, resulting in the small SHG signal from the

bilayer (see the main text for discussion). Accounting for the electronic coupling between the layers may result in additional contributions to the SHG which has the same smallness $\propto q_\omega a$.

In Eqs. (4) we neglected the differences between the transmission coefficients of the monolayers and unity at both fundamental and double frequencies and disregarded the effect of the substrate on the fields. While the transmission coefficient of the monolayer $t_\omega$ at a fundamental frequency is very close to 1 because the incident radiation is non-resonant, the effect of the monolayers on the transmission and reflection of the second harmonic can be sizeable due to the exciton resonances. With account for multiple reflections of the second harmonics from the monolayers, we have an enhancement of both the reflected and transmitted waves. In the reflection geometry the enhancement factor [cf. Ref. [66]]

$$F_r = 1 + t_{2\omega} \frac{r_{2\omega} e^{i2q_{2\omega}a}}{1 - r_{2\omega}^2 e^{i2q_{2\omega}a}} \approx \frac{1}{1 - r_{2\omega}}, \quad (5)$$

where $r_{2\omega}$ is the reflection coefficient from the monolayer at the frequency $2\omega$, $t_{2\omega} = 1 + r_{2\omega}$, and in the last equation we neglected the phase acquired by the field propagating between the layers. Making use of the explicit form of the monolayer reflection coefficient we recast $F_r$ as

$$F_r = \frac{\omega_0 - 2\omega - i(\Gamma_0 + \Gamma)}{\omega_0 - 2\omega - i(2\Gamma_0 + \Gamma)},$$

where $\omega_0$ is the exciton resonance frequency, $\Gamma_0$ and $\Gamma$ are its radiative and non-radiative damping.

The expression for the enhancement factor in the transmission geometry has somewhat more complex form

$$F_t = \frac{e^{iq_{2\omega}a} \frac{t_{2\omega}}{1-r_{2\omega}^2 e^{iq_{2\omega}a}} - e^{i2q_{\omega}a} \left(1 + \frac{t_{2\omega}r_{2\omega}e^{2iq_{2\omega}a}}{1-r_{2\omega}^2 e^{2iq_{2\omega}a}}\right)}{t_{2\omega}e^{iq_{2\omega}a} - e^{i2q_{\omega}a}}. \quad (6)$$

It follows from Eqs. (5) and (6) that in the vicinity of the exciton resonance the enhancement due to multiple reflections is unimportant.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 8: Comparison between SHG and differential reflectivity spectra of monolayer and trilayer MoS$_2$.</figcaption></figure>

It is instructive to discuss the relation between the effects described above and the surface SHG effect both at oblique and normal incidence studied in a number of works, including Refs. [52, 55, 56, 67, 68]. It has been previously established that either the surface of a media has a reduced symmetry, e.g., [56] and allows for the SHG by itself, or the SHG arises due to the quadrupolar or magneto-dipolar transitions ($qE^2$ or EB mechanisms), e.g., [55]. The mechanism discussed above belongs to the second type of the effects.

To conclude this part, in ideal centrosymmetric bilayer the SHG is possible due to the effects related to the finite wavevector of the radiation and a finite width of

<figure><img src="image_13.png" /><figcaption>FIG. 9: Comparison between SHG of bilayer MoS$_2$ with hBN. Notice that the SHG from hBN is multiplied by a factor of 50.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 10: Power dependent SHG in resonance with the A-exciton clearly demonstrates a two-photon process.</figcaption></figure>

the bilayer. Because of this, the effect is suppressed in pristine bilayer as compared to the non-centrosymmetric monolayer.

### Electric field induced mechanism

#### Minimal model

Let us first consider the minimal model demonstrating the electric field induced SHG in bilayers. To that end motivated by microscopic analysis presented in Refs. [40, 47] we consider the model of two interlayer exciton states $IE_{1,2}$ and two intralayer exciton states $B_{1,2}$. These four
states are close in energy and are effectively mixed by the hole tunneling. Let $T$ be the (hole) tunneling constant which mixes corresponding IE and B exciton states. The Hamiltonian of the system in the presence of the electric field $F_z$ reads (the basic states are $IE_1$, $B_1$, $IE_2$, $B_2$):

$$ \mathcal{H} = \begin{pmatrix} E_{\text{IE}} - \mu F_z & T & 0 & 0 \\ T^* & E_B & 0 & 0 \\ 0 & 0 & E_{\text{IE}} + \mu F_z & T \\ 0 & 0 & T^* & E_B \end{pmatrix}. \quad (7) $$

Here $\mu$ is the indirect exciton static dipole matrix element along the structure normal (see main text for details).

For each pair of levels $IE_1$ & $B_1$ and $IE_2$ & $B_2$ we have a $2 \times 2$ Hamiltonian which is readily diagonalized as

$$ \begin{aligned} E_{1,\pm}(F_z) &= \frac{E_{\text{IE}} - \mu F_z + E_B}{2} & (8a) \\ &\pm \sqrt{\left(\frac{E_{\text{IE}} - \mu F_z - E_B}{2}\right)^2 + |T|^2}, \\ E_{2,\pm}(F_z) &= E_{1,\pm}(-F_z), & (8b) \end{aligned} $$

with the eigenfunctions

$$ |1, \pm\rangle = \alpha_{\pm}(F_z)|\text{IE}_1\rangle + \beta_{\pm}(F_z)|B_1\rangle, \quad (8c) $$

$$ |2, \pm\rangle = \alpha_{\pm}(-F_z)|\text{IE}_2\rangle + \beta_{\pm}(-F_z)|B_2\rangle. \quad (8d) $$

The coefficients $\alpha_{\pm}(F_z)$, $\beta_{\pm}(F_z)$ read

$$ \begin{aligned} \alpha_{\pm}(F_z) &= \frac{E_{1,\pm} - E_B}{T^*} \beta_{\pm}(F_z), & (8e) \\ \beta_{\pm}(F_z) &= \frac{1}{\sqrt{1 + \frac{(E_{1,\pm} - E_B)^2}{|T|^2}}}. \end{aligned} $$

In order to derive compact analytical expressions we assume that the distance between the IE and B excitons $\Delta = E_B - E_{\text{IE}}$ exceeds by far both $|T|$ and $|\mu F_z|$, and also that $|\mu F_z| \ll |T|$. Hence, the eigenenergies and eigenstates are

$$ \widetilde{\text{IE}_1}: \quad E_{\widetilde{\text{IE}_1}} \approx E_{\text{IE}} - \mu F_z - \frac{|T|^2}{\Delta}, \quad |\widetilde{\text{IE}_1}\rangle = |\text{IE}_1\rangle - \frac{T^*}{\Delta} \left(1 - \frac{\mu F_z}{\Delta}\right) |B_1\rangle, \quad (9a) $$

$$ \widetilde{\text{IE}_2}: \quad E_{\widetilde{\text{IE}_2}} \approx E_{\text{IE}} + \mu F_z - \frac{|T|^2}{\Delta}, \quad |\widetilde{\text{IE}_2}\rangle = |\text{IE}_2\rangle - \frac{T^*}{\Delta} \left(1 + \frac{\mu F_z}{\Delta}\right) |B_2\rangle, \quad (9b) $$

$$ \widetilde{\text{B}_1}: \quad E_{\widetilde{\text{B}_1}} \approx E_B + \frac{|T|^2}{\Delta} \left(1 - \frac{\mu F_z}{\Delta}\right), \quad (9c) \\ |\widetilde{\text{B}_1}\rangle = \left[1 - \frac{|T|^2}{2\Delta^2} \left(1 - \frac{2\mu F_z}{\Delta}\right)\right] |\text{B}_1\rangle + \frac{T}{\Delta} \left(1 - \frac{\mu F_z}{\Delta}\right) |\text{IE}_1\rangle, $$

$$ \widetilde{\text{B}_2}: \quad E_{\widetilde{\text{B}_2}} \approx E_B + \frac{|T|^2}{\Delta} \left(1 + \frac{\mu F_z}{\Delta}\right), \quad (9d) $$

$$ |\widetilde{\text{B}_2}\rangle = \left[1 - \frac{|T|^2}{2\Delta^2} \left(1 + \frac{2\mu F_z}{\Delta}\right)\right] |\text{B}_2\rangle + \frac{T}{\Delta} \left(1 + \frac{\mu F_z}{\Delta}\right) |\text{IE}_2\rangle. $$

Let $M_2$ and $M_1$ be the two-photon and the one-photon excitation matrix elements of the B-exciton. The susceptibilities of $B_{1,2}$ excitons in isolated monolayers $T=0$ read ($C$ is a constant) [36]

$$ \chi_{\text{B},1,2} = \pm C \frac{M_1 M_2}{2\hbar\omega - E_B + i\Gamma_B}. \quad (10) $$

For a bilayer at $F_z=0$ the sum $\chi_{\text{B},1} + \chi_{\text{B},2} = 0$ signifying the presence of an inversion center.

The SHG generation results from the mixing of those states with B-excitons. Generally, the susceptibility can be written as

$$ \chi = -CM_1 M_2 \sum_{j=1}^{2} \sum_{\pm} (-1)^j \frac{|\beta_{\pm}(F_z)|^2}{2\hbar\omega - E_{j,\pm}(F_z) + i\Gamma_{j,\pm}}. \quad (11) $$

As a result, for the weak tunneling and field we obtain from Eqs. (9)
$$ \widetilde{\text{IE}_1}: \quad \chi_{\text{IE}_1} = \frac{|T|^2}{\Delta^2} \left(1 - \frac{2\mu F_z}{\Delta}\right) \frac{CM_1 M_2}{2\hbar\omega - E_{\text{IE}} + \mu F_z + \frac{|T|^2}{\Delta} + i\Gamma_{\text{IE}}}, \quad (12a) $$

$$ \widetilde{\text{IE}_2}: \quad \chi_{\text{IE}_2} = -\frac{|T|^2}{\Delta^2} \left(1 + \frac{2\mu F_z}{\Delta}\right) \frac{CM_1 M_2}{2\hbar\omega - E_{\text{IE}} - \mu F_z + \frac{|T|^2}{\Delta} + i\Gamma_{\text{IE}}}, \quad (12b) $$

$$ \widetilde{\text{B}_1}: \quad \chi_{\text{B}_1} = \left[1 - \frac{|T|^2}{2\Delta^2} \left(1 - \frac{2\mu F_z}{\Delta}\right)\right] \frac{CM_1 M_2}{2\hbar\omega - E_{\text{B}} - \frac{|T|^2}{\Delta} \left(1 - \frac{\mu F_z}{\Delta}\right) + i\Gamma_{\text{B}}}, \quad (12c) $$

$$ \widetilde{\text{B}_2}: \quad \chi_{\text{B}_2} = -\left[1 - \frac{|T|^2}{2\Delta^2} \left(1 + \frac{2\mu F_z}{\Delta}\right)\right] \frac{CM_1 M_2}{2\hbar\omega - E_{\text{B}} - \frac{|T|^2}{\Delta} \left(1 + \frac{\mu F_z}{\Delta}\right) + i\Gamma_{\text{B}}}. \quad (12d) $$

The total susceptibility

$$ \chi = \chi_{\text{IE}_1} + \chi_{\text{IE}_2} + \chi_{\text{B}_1} + \chi_{\text{B}_2}. \quad (13) $$

It vanishes at $F_z = 0$. At $F_z \neq 0$ the inversion symmetry is broken and effect arises both due to the energy shifts of the excitons and due to the linear in $F_z$ contributions to the oscillator strength.

To get better insight into the SHG in bilayer materials and for the comparison with available experiments it is instructive to analyze the situation where the electric field is so weak such that $|\mu F_z| \ll \Gamma_{\text{B,IE},...}$. The susceptibility in the vicinity of the IE excitons takes the form

$$ \begin{aligned} \chi_{\text{IE}} &= \chi_{\text{IE}_1} + \chi_{\text{IE}_2} \\ &= -\frac{2\mu F_z}{2\hbar\omega - E_{\text{IE}} + \frac{|T|^2}{\Delta} + i\Gamma_{\text{IE}}} \frac{CM_1 M_2 |T|^2}{\Delta^2} \times \\ & \quad \left(\frac{2}{\Delta} + \frac{1}{2\hbar\omega - E_{\text{IE}} + \frac{|T|^2}{\Delta} + i\Gamma_{\text{IE}}}\right). \end{aligned} \quad (14) $$

For reasonable $\Delta > \Gamma_{\text{IE}}$ the second term in parentheses dominates in the vicinity of the resonance. It results in Eq. (2) of the main text with the notation $C_2|\Phi_{\text{B:1s}}(0)|^2 = CM_1 M_2$, cf. Ref. [36]. Similar expression can be derived describing the susceptibility in the vicinity of the B excitons:

$$ \begin{aligned} \chi_B &= \chi_{\text{B}_1} + \chi_{\text{B}_2} \\ &= -\frac{2\mu F_z}{2\hbar\omega - E_B - \frac{|T|^2}{\Delta} + i\Gamma_B} \frac{CM_1 M_2 |T|^2}{\Delta^2} \times \\ & \quad \left(\frac{1}{\Delta} - \frac{1}{2\hbar\omega - E_B - \frac{|T|^2}{\Delta} + i\Gamma_B}\right). \end{aligned} \quad (15) $$

### Extensions of the model

The deficiency of the model above is the absence of the effect on A-excitons and also the neglect of the mixing of

intralayer excitons due to the tunneling and symmetry-breaking perturbations.

Let us first include exciton tunneling as a whole described by the matrix element $T_B$ in the Hamiltonian, which tends to form symmetric and antisymmetric combinations of B-excitons as well as symmetry breaking between the MLs in bilayer resulting in the splitting of B-excitons with $E_{B,1} - E_{B,2} = \delta E_B$:

$$ \mathcal{H} = \begin{pmatrix} E_{\text{IE}} - \mu F_z & T & 0 & 0 \\ T^* & E_{B,1} & 0 & T_B \\ 0 & 0 & E_{\text{IE}} + \mu F_z & T \\ 0 & T_B^* & T^* & E_{B,2} \end{pmatrix}. \quad (16) $$

We focus at the reasonable relation between the parameters of the system with the electric field being low as compared with the broadening of the indirect exciton, $|\mu F_z| \ll \Gamma_{\text{IE}} \ll \Delta$, but, on the other hand, sufficiently high so that $|\mu F_z| \gg |TT_B/\Delta|$, $\delta E_B |T|^2/\Delta^2$. Under these approximations the $T_B$ and $\delta E_B$ do not affect the IE exciton and its susceptibility is given by Eq. (14).

The effective Hamiltonian describing the doublet of B-excitons in the second-order in $T$ and lowest order in $\mu F_z$, $\delta E_B$ reads

$$ \mathcal{H} = \begin{pmatrix} E_{B,1} + \mu F_z \frac{|T|^2}{\Delta^2} + \frac{|T|^2}{\Delta} & T_B \\ T_B^* & E_{B,2} - \mu F_z \frac{|T|^2}{\Delta^2} + \frac{|T|^2}{\Delta} \end{pmatrix}. \quad (17) $$

Furthermore, in a weak asymmetry limit where $|T_B| \gg \delta E_B$, $|\mu F_z||T|^2/\Delta^2$, the Hamiltonian (17) in the basis of $(|B_1\rangle \pm |B_2\rangle)/\sqrt{2}$ states transforms to (we assume that the diagonal terms $|T|^2/\Delta$ are included in $E_B$)

$$ \mathcal{H} = \begin{pmatrix} E_B + T_B & \delta E_B + \mu F_z \frac{|T|^2}{\Delta^2} \\ \delta E_B + \mu F_z \frac{|T|^2}{\Delta^2} & E_B - T_B \end{pmatrix}, \quad (18) $$
whose eigenstates are:

$$|B_+\rangle: \quad E_{B+} = E_B + T_B, \quad (19)$$

$$|B_+\rangle = \frac{|B_1\rangle + |B_2\rangle}{\sqrt{2}} + \frac{\delta}{2T_B} \frac{|B_1\rangle - |B_2\rangle}{\sqrt{2}},$$

$$|B_-\rangle: \quad E_{B-} = E_B - T_B, \quad (20)$$

$$|B_+\rangle = \frac{|B_1\rangle - |B_2\rangle}{\sqrt{2}} - \frac{\delta}{2T_B} \frac{|B_1\rangle + |B_2\rangle}{\sqrt{2}},$$

with $\delta = \delta E_B + \mu F_z |T|^2/\Delta^2$. As a result, the non-linear susceptibility in the spectral range of B-excitons reads

$$\chi = \frac{\delta}{2T_B} C M_1 M_2 \times \left( \frac{1}{2\hbar\omega - E_B - T_B + i\Gamma_B} - \frac{1}{2\hbar\omega - E_B + T_B + i\Gamma_B} \right) \approx \left( \delta E_B + \frac{\mu F_z |T|^2}{\Delta^2} \right) \frac{C M_1 M_2}{(2\hbar\omega - E_B + i\Gamma_B)^2 - T_B^2} \quad (21)$$

Very similar treatment can be made for A-exciton assuming that it is mixed with the IE(B) exciton, which follows from the spin conservation rule. Making natural changes: $T_B$ to the tunneling matrix element of the A-exciton $T_A$, $\Delta$ to the splitting $\Delta' = E_A - E_{IE(B)}$, and $\delta E_B$ to the A-exciton splitting $E_{A,1} - E_{A,2} = \delta E_A$, we obtain

$$\chi_A = \delta' \frac{C' M_1 M_2}{(2\hbar\omega - E_A + i\Gamma_A)^2 - T_A^2}. \quad (22)$$

Here $\delta' = \delta E_A + \mu F_z |T_A|^2/\Delta'^2$. It follows from Eq. (22) that the effect of the electric field on A-exciton is weaker than on the IE. It follows from atomistic calculations that the A-excitons can be mixed with interlayer exciton $IE_1$ via the spin-flip electron tunneling. This effect can be included in our model in a same way, resulting in susceptibilities similar to Eqs. (14) and (15) derived for the mixing of the IE- and B-excitons.

**Acknowledgements.**— Toulouse acknowledges funding from ANR 2D-vdW-Spin, ANR MagicValley, ANR IXTASE, ANR HiLight, ITN 4PHOTON Marie Sklodowska Curie Grant Agreement No. 721394 and the Institut Universitaire de France. Growth of hexagonal boron nitride crystals was supported by the Elemental Strategy Initiative conducted by the MEXT, Japan, Grant Number JPMXP0112101001, JSPS KAKENHI Grant Number JP20H00354 and the CREST(JPMJCR15F3), JST. M.M.G. and L.E.G. acknowledge the RFBR and CNRS joint project 20-52-16303. L.E.G. was supported by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”.

* Electronic address: glazov@coherent.ioffe.ru

† Electronic address: urbaszek@insa-toulouse.fr

[1] Zhang, J., Zhao, W., Yu, P., Yang, G. & Liu, Z. Second harmonic generation in 2D layered materials. 2D Materials 7, 042002 (2020).

[2] Trovatello, C. et al. Optical parametric amplification by monolayer transition metal dichalcogenides. Nature Photonics 15, 6–10 (2020).

[3] Franken, P. A., Hill, A. E., Peters, C. W. & Weinreich, G. Generation of optical harmonics. Phys. Rev. Lett. 7, 118–119 (1961).

[4] Terhune, R. W., Maker, P. D. & Savage, C. M. Optical harmonic generation in calcite. Phys. Rev. Lett. 8, 404–406 (1962).

[5] Yariv, A. & Yeh, P. Photonics: optical electronics in modern communications, vol. 6 (Oxford university press New York, 2007).

[6] Heinz, T. F., Chen, C., Ricard, D. & Shen, Y. Spectroscopy of molecular monolayers by resonant second-harmonic generation. Physical Review Letters 48, 478 (1982).

[7] Shree, S., Paradisanos, I., Marie, X., Robert, C. & Urbaszek, B. Guide to optical spectroscopy of layered semiconductors. Nature Reviews Physics 3, 39–54 (2021).

[8] Fiebig, M., Pavlov, V. V. & Pisarev, R. V. Second-harmonic generation as a tool for studying electronic and magnetic structures of crystals. JOSA B 22, 96–118 (2005).

[9] Lafrentz, M. et al. Second-harmonic generation spectroscopy of excitons in ZnO. Physical Review B 88, 235207 (2013).

[10] Bergfeld, S. & Daum, W. Second-harmonic generation in GaAs: Experiment versus theoretical predictions of $\chi_{xyz}^{(2)}$. Physical review letters 90, 036801 (2003).

[11] Yakovlev, D. R. et al. Exciton spectroscopy of semiconductors by the method of optical harmonics generation (review). Physics of the Solid State 60, 1471–1486 (2018).

[12] De Dominicis, L. et al. Second- and third- harmonic generation in single-walled carbon nanotubes at nanosecond time scale. Applied Physics Letters 85, 1418–1420 (2004).

[13] Sun, Z. et al. Giant nonreciprocal second-harmonic generation from antiferromagnetic bilayer CrI$_3$. Nature 572, 497–501 (2019).

[14] Löchner, F. J. et al. Controlling second-harmonic diffraction by nano-patterning MoS$_2$ monolayers. Optics Express 27, 35475–35484 (2019).

[15] Cunha, R. et al. Second harmonic generation in defective hexagonal boron nitride. Journal of Physics: Condensed Matter 32, 19LT01 (2020).

[16] Seyler, K. L. et al. Electrical control of second-harmonic generation in a WSe$_2$ monolayer transistor. Nature nanotechnology 10, 407–411 (2015).

[17] Wang, G. et al. Giant enhancement of the optical second-harmonic emission of WSe$_2$ monolayers by laser excitation at exciton resonances. Physical review letters 114, 097403 (2015).

[18] Lin, K.-Q., Bange, S. & Lupton, J. M. Quantum interference in second-harmonic generation from monolayer WSe$_2$. Nature Physics 15, 242–246 (2019).

[19] Mennel, L. et al. Optical imaging of strain in two-
dimensional crystals. *Nature communications* **9**, 1–6 (2018).

[20] Abdelwahab, I. et al. Highly enhanced third-harmonic generation in 2D perovskites at excitonic resonances. *ACS nano* **12**, 644–650 (2018).

[21] Fiebig, M., Frohlich, D., Leute, S. & Pisarev, R. V. Second harmonic spectroscopy and control of domain size in antiferromagnetic YMnO₃. *Journal of Applied Physics* **83**, 6560–6562 (1998).

[22] Glazov, M. M. & Ganichev, S. D. High frequency electric field induced nonlinear effects in graphene. *Physics Reports* **535**, 101–138 (2014).

[23] Dean, J. J. & van Driel, H. M. Second harmonic generation from graphene and graphitic films. *Applied Physics Letters* **95**, 261910 (2009).

[24] Hendry, E., Hale, P. J., Moger, J., Savchenko, A. & Mikhailov, S. A. Coherent nonlinear optical response of graphene. *Physical review letters* **105**, 097401 (2010).

[25] Mikhailov, S. A. Theory of the giant plasmon-enhanced second-harmonic generation in graphene and semiconductor two-dimensional electron systems. *Phys. Rev. B* **84**, 045432 (2011).

[26] Glazov, M. M. Second harmonic generation in graphene. *JETP letters* **93**, 366–371 (2011).

[27] Gullans, M., Chang, D., Koppens, F., de Abajo, F. G. & Lukin, M. D. Single-photon nonlinear optics with graphene plasmons. *Physical review letters* **111**, 247401 (2013).

[28] Soavi, G. et al. Hot electrons modulation of third-harmonic generation in graphene. *ACS Photonics* **6**, 2841–2849 (2019).

[29] Yu, R., Cox, J. D. & de Abajo, F. J. G. Nonlinear plasmonic sensing with nanographene. *Phys. Rev. Lett.* **117**, 123904 (2016).

[30] Golub, L. & Tarasenko, S. Valley polarization induced second harmonic generation in graphene. *Physical Review B* **90**, 201402 (2014).

[31] Bykov, A. Y., Murzina, T. V., Rybin, M. G. & Obraztsova, E. D. Second harmonic generation in multilayer graphene induced by direct electric current. *Phys. Rev. B* **85**, 121413 (2012).

[32] Li, Y. et al. Probing symmetry properties of few-layer MoS₂ and h-BN by optical second-harmonic generation. *Nano letters* **13**, 3329–3333 (2013).

[33] Malard, L. M., Alencar, T. V., Barboza, A. P. M., Mak, K. F. & de Paula, A. M. Observation of intense second harmonic generation from MoS₂ atomic crystals. *Phys. Rev. B* **87**, 201401 (2013).

[34] Wang, G. et al. Colloquium: Excitons in atomically thin transition metal dichalcogenides. *Reviews of Modern Physics* **90**, 021001 (2018).

[35] Wang, G. et al. Exciton states in monolayer MoSe₂: impact on interband transitions. *2D Materials* **2**, 045005 (2015).

[36] Glazov, M. et al. Intrinsic exciton-state mixing and nonlinear optical properties in transition metal dichalcogenide monolayers. *Physical Review B* **95**, 035311 (2017).

[37] Zhao, M. et al. Atomically phase-matched second-harmonic generation in a 2D crystal. *Light: Science & Applications* **5**, e16131–e16131 (2016).

[38] Wen, X., Gong, Z. & Li, D. Nonlinear optics of two-dimensional transition metal dichalcogenides. *InfoMat* **1**, 317–337 (2019).

[39] Weismann, M. & Panoiu, N. C. Theoretical and computational analysis of second-and third-harmonic generation in periodically patterned graphene and transition-metal dichalcogenide monolayers. *Physical Review B* **94**, 035435 (2016).

[40] Leisgang, N., Shree, S., Paradisanos, I. et al. Giant stark splitting of an exciton in bilayer MoS₂. *Nat. Nanotechnol.* **15**, 901–907 (2020).

[41] Klein, J. et al. Electric-field switchable second-harmonic generation in bilayer MoS₂ by inversion symmetry breaking. *Nano Letters* **17**, 392–398 (2017).

[42] Yu, H., Talukdar, D., Xu, W., Khurgin, J. B. & Xiong, Q. Charge-induced second-harmonic generation in bilayer WSe₂. *Nano letters* **15**, 5653–5657 (2015).

[43] Deilmann, T. & Thygesen, K. S. Interlayer excitons with large optical amplitudes in layered van der Waals materials. *Nano letters* **18**, 2984–2989 (2018).

[44] Slobodeniuk, A. et al. Fine structure of k-excitons in multilayers of transition metal dichalcogenides. *2D Materials* **6**, 025026 (2019).

[45] Gerber, I. C. et al. Interlayer excitons in bilayer MoS₂ with strong oscillator strength up to room temperature. *Physical Review B* **99**, 035443 (2019).

[46] Paradisanos, I. et al. Controlling interlayer excitons in MoS₂ layers grown by chemical vapor deposition. *Nature communications* **11**, 1–7 (2020).

[47] Lorchat, E. et al. Excitons in bilayer MoS₂ displaying a colossal electric field splitting and tunable magnetic response. *Phys. Rev. Lett.* **126**, 037401 (2021).

[48] Taniguchi, T. & Watanabe, K. Synthesis of high-purity boron nitride single crystals under high pressure by using Ba–BN solvent. *Journal of crystal growth* **303**, 525–529 (2007).

[49] Stier, A. V., Wilson, N. P., Clark, G., Xu, X. & Crooker, S. A. Probing the influence of dielectric environment on excitons in monolayer WSe₂: insight from high magnetic fields. *Nano letters* **16**, 7054–7060 (2016).

[50] Raja, A. et al. Dielectric disorder in two-dimensional materials. *Nature nanotechnology* **14**, 832–837 (2019).

[51] Rhodes, D., Chae, S. H., Ribeiro-Palau, R. & Hone, J. Disorder in van der Waals heterostructures of 2D materials. *Nature materials* **18**, 541 (2019).

[52] Shen, Y. Surface properties probed by second-harmonic and sum-frequency generation. *Nature* **337**, 519–525 (1989).

[53] Cadiz, F. et al. Excitonic linewidth approaching the homogeneous limit in MoS₂-based van der Waals heterostructures. *Physical Review X* **7**, 021026 (2017).

[54] Goryca, M. et al. Revealing exciton masses and dielectric properties of monolayer semiconductors with high magnetic fields. *Nature communications* **10**, 1–12 (2019).

[55] Bloembergen, N., Chang, R. K., Jha, S. S. & Lee, C. H. Optical second-harmonic generation in reflection from media with inversion symmetry. *Phys. Rev.* **174**, 813–822 (1968).

[56] Maki, J. J., Kauranen, M. & Persoons, A. Surface second-harmonic generation from chiral materials. *Phys. Rev. B* **51**, 1425–1434 (1995).

[57] Stiehm, T. et al. Supercontinuum second harmonic generation spectroscopy of atomically thin semiconductors. *Review of Scientific Instruments* **90**, 083102 (2019).

[58] Gong, Z. et al. Magnetoelectric effects and valley-controlled spin quantum gates in transition metal dichalcogenide bilayers. *Nature Communications* **4**, 2053 (2013).
[59] Roch, J. G. et al. Quantum-confined stark effect in a MoS₂ monolayer van der Waals heterostructure. *Nano letters* **18**, 1070–1074 (2018).

[60] Verzhbitskiy, I., Vella, D., Watanabe, K., Taniguchi, T. & Eda, G. Suppressed out-of-plane polarizability of free excitons in monolayer WSe₂. *ACS nano* **13**, 3218–3224 (2019).

[61] Horng, J. et al. Observation of interlayer excitons in mose₂ single crystals. *Phys. Rev. B* **97**, 241404 (2018).

[62] Alexeev, E. M. et al. Resonantly hybridized excitons in moiré superlattices in van der waals heterostructures. *Nature* **567**, 81–86 (2019).

[63] Tang, Y. et al. Tuning layer-hybridized moiré excitons by the quantum-confined stark effect. *Nature Nanotechnology* **16**, 52–57 (2021).

[64] Refaely-Abramson, S., Qiu, D. Y., Louie, S. G. & Neaton, J. B. Defect-induced modification of low-lying excitons and valley selectivity in monolayer transition metal dichalcogenides. *Physical review letters* **121**, 167402 (2018).

[65] Robert, C. et al. Optical spectroscopy of excited exciton states in mos 2 monolayers in van der waals heterostructures. *Physical Review Materials* **2**, 011001 (2018).

[66] Fang, H. H. et al. Control of the exciton radiative lifetime in van der waals heterostructures. *Phys. Rev. Lett.* **123**, 067401 (2019).

[67] Guyot-Sionnest, P., Chen, W. & Shen, Y. R. General considerations on optical second-harmonic generation from surfaces and interfaces. *Phys. Rev. B* **33**, 8254–8263 (1986).

[68] Sipe, J. E., Moss, D. J. & van Driel, H. M. Phenomenological theory of optical second- and third-harmonic generation from cubic centrosymmetric crystals. *Phys. Rev. B* **35**, 1129–1141 (1987).