Non-adiabatic exciton-phonon coupling in Raman spectroscopy of layered materials

TL;DR

This paper introduces an ab initio computational method to analyze resonant Raman intensities in layered materials, accounting for excitonic and non-adiabatic effects, revealing quantum interference phenomena affecting exciton behavior.

Contribution

It develops a novel diagrammatic approach for calculating resonant Raman intensities including excitonic and non-adiabatic effects in layered materials.

Findings

Quantum interference causes strong redistribution of oscillator strength in h-BN.

Spin-orbit splitting suppresses quantum interference effects in MoS2.

The approach enables detailed analysis of exciton-phonon interactions in layered materials.

Abstract

We present an ab initio computational approach for the calculation of resonant Raman intensities, including both excitonic and non-adiabatic effects. Our diagrammatic approach, which we apply to two prototype, semiconducting layered materials, allows a detailed analysis of the impact of phonon-mediated exciton-exciton scattering on the intensities. In the case of bulk hexagonal boron nitride, this scattering leads to strong quantum interference between different excitonic resonances, strongly redistributing oscillator strength with respect to optical absorption spectra. In the case of MoS$_2$, we observe that quantum interference effects are suppressed by the spin-orbit splitting of the excitons.

Equations

$${\dot{P}}_{1-ph.}=\text{const.}\times {\omega }_{\mathrm{L}}{\omega }_{\mathrm{D}}\sum _{\lambda }{\left|{\tilde{\mathcal{M}}}_{\mu \nu }^{\lambda }({\omega }_{\mathrm{L}}-{\omega }_{\mathrm{D}},{\omega }_{\mathrm{D}})\right|}^2{\mathcal{A}}_{\lambda }({\omega }_{\mathrm{L}}-{\omega }_{\mathrm{D}}).$$

$$\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}(\omega^{\prime},\omega)\equiv\frac{1}{\hbar^{2}}\int_{-\infty}^{+\infty}\mathrm{d}t\,\mathrm{e}^{i\omega t}\int_{-\infty}^{+\infty}\mathrm{d}t^{\prime}\,\mathrm{e}^{i\omega^{\prime}t^{\prime}}(-i)\langle 0|\mathcal{T}\left[\hat{F}_{\lambda}(t^{\prime})\hat{J}_{\nu}(t)\hat{J}^{\dagger}_{\mu}(0)\right]|0\rangle\big{|}_{\mathrm{connect.}},$$

$$\begin{split}\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}(\omega_{\lambda},\omega_{\mathrm{L}}-\omega_{\lambda})&=\sum_{S,S^{\prime}}\frac{d^{\mu}_{S}\big{(}g^{\lambda}_{S,S^{\prime}}\big{)}^{*}\big{(}d^{\nu}_{S^{\prime}}\big{)}^{*}}{(\hbar\omega_{\mathrm{L}}-E_{S}+i\gamma)(\hbar\omega_{\mathrm{L}}-\hbar\omega_{\lambda}-E_{S^{\prime}}+i\gamma)}\\ &\quad+\sum_{S,S^{\prime}}\frac{\big{(}d^{\mu}_{S}\big{)}^{*}g^{\lambda}_{S,S^{\prime}}d^{\nu}_{S^{\prime}}}{(\hbar\omega_{\mathrm{L}}+E_{S}-i\gamma)(\hbar\omega_{\mathrm{L}}-\hbar\omega_{\lambda}+E_{S^{\prime}}-i\gamma)}.\end{split}$$

$$d^{\mu}_{S}\equiv\sum_{\mathbf{k},\begin{subarray}{c}a\in\mathcal{C}\\ b\in\mathcal{V}\end{subarray}}d^{\mu,(\mathrm{b})}_{\mathbf{k},a,b}\big{(}A^{S}_{\mathbf{k},a,b}\big{)}^{*}$$

$$g^{\lambda}_{S^{\prime},S}\equiv\sum_{\mathbf{k}}\bigg{[}\sum_{\begin{subarray}{c}a,c\in\mathcal{C}\\ b\in\mathcal{V}\end{subarray}}\big{(}A^{S^{\prime}}_{\mathbf{k},a,b}\big{)}^{*}g^{\lambda}_{\mathbf{k},a,c}A^{S}_{\mathbf{k},c,b}-\sum_{\begin{subarray}{c}a\in\mathcal{C}\\ b,c\in\mathcal{V}\end{subarray}}\big{(}A^{S^{\prime}}_{\mathbf{k},a,b}\big{)}^{*}g^{\lambda}_{\mathbf{k},c,b}A^{S}_{\mathbf{k},a,c}\bigg{]},$$

$$\begin{split}\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}\big{|}_{\text{intra}}(\omega_{\lambda},\omega_{\mathrm{L}}-\omega_{\lambda})=\sum_{S,S^{\prime}}\bigg{|}_{E_{S}=E_{S^{\prime}}}\bigg{[}&\frac{d^{\mu}_{S}\big{(}g^{\lambda}_{S,S^{\prime}}\big{)}^{*}\big{(}d^{\nu}_{S^{\prime}}\big{)}^{*}}{(\hbar\omega_{\mathrm{L}}-E_{S}+i\gamma)(\hbar\omega_{\mathrm{L}}-\hbar\omega_{\lambda}-E_{S^{\prime}}+i\gamma)}\\ &+\frac{\big{(}d^{\mu}_{S}\big{)}^{*}g^{\lambda}_{S,S^{\prime}}d^{\nu}_{S^{\prime}}}{(\hbar\omega_{\mathrm{L}}+E_{S}-i\gamma)(\hbar\omega_{\mathrm{L}}-\hbar\omega_{\lambda}+E_{S^{\prime}}-i\gamma)}\bigg{]}.\end{split}$$

$$\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}\big{|}_{\text{inter}}(\omega_{\lambda},\omega_{\mathrm{L}}-\omega_{\lambda})\equiv\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}(\omega_{\lambda},\omega_{\mathrm{L}}-\omega_{\lambda})-\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}\big{|}_{\text{intra}}(\omega_{\lambda},\omega_{\mathrm{L}}-\omega_{\lambda})$$

Full text

Non-adiabatic exciton-phonon coupling in Raman spectroscopy of layered materials

Sven Reichardt,1∗ Ludger Wirtz1

1Department of Physics and Materials Science, University of Luxembourg

1511 Luxembourg, Luxembourg

∗To whom correspondence should be addressed; E-mail: sven.reichardt@uni.lu

We present an ab initio computational approach for the calculation of resonant Raman intensities, including both excitonic and non-adiabatic effects. Our diagrammatic approach, which we apply to two prototype, semiconducting layered materials, allows a detailed analysis of the impact of phonon-mediated exciton-exciton scattering on the intensities. In the case of bulk hexagonal boron nitride, this scattering leads to strong quantum interference between different excitonic resonances, strongly redistributing oscillator strength with respect to optical absorption spectra. In the case of MoS2, we observe that quantum interference effects are suppressed by the spin-orbit splitting of the excitons.

Introduction

Raman spectroscopy is a versatile tool for the characterization of many materials and, in particular, 2D materials  (?). It allows the study of both vibrational and electronic properties. Among many other things, it serves to probe many-body  (?) and non-adiabatic effects  (?, ?) as well as excitonic ones  (?, ?). From a theoretical point of view, however, a general and comprehensive approach that captures both of these effects is missing so far. Regarding one-phonon Raman scattering, recent theoretical and computational efforts have focused either on a leading-order perturbation theory approach   (?, ?, ?, ?, ?, ?, ?), which is able to capture non-adiabatic effects but not excitonic ones, or on a method based on the calculation of the change of the dielectric susceptibility with static atomic displacements along phonon modes  (?, ?, ?, ?, ?, ?), which captures excitonic effects but is entirely adiabatic. As such, an approach that is able to capture both effects is highly desirable not only from a theoretical and conceptional point of view, but also in practice. This is especially true for systems with sizable phonon frequencies and strong excitonic effects or also for materials whose properties are a priori unknown.

In a recent work  (?), we laid the groundwork for such a comprehensive theoretical description of Raman scattering from first principles that allows the inclusive description of both excitonic and non-adiabatic effects in one fully quantum mechanical theory. Here we now apply this general theory to the case of resonant, one-phonon Raman scattering. We present concrete expressions for the Raman intensity in terms of quantities that can be obtained by ab initio calculations on the level of density functional theory (DFT) and many-body perturbation theory (MBPT), i.e., within the $GW$-approximation for the quasi-particle band structure and solving the Bethe-Salpeter equation to include excitonic effects. We then apply it to two materials, bulk hexagonal boron nitride (hBN) and monolayer molybdenum disulfide (MoS2), which are known to feature strong excitonic effects in both absorption  (?, ?) and Raman spectra  (?, ?) due to their layered, quasi-two-dimensional geometry. The modular structure of our method allows us to further analyze our numerical results in detail and explain them in a physically intuitive way.

In the case of hBN our two main findings are (i) the suppression of some of the higher excitonic resonances in the Raman spectrum, contrary to their brightness in the optical absorption spectrum, and (ii) strong quantum interference between the first two excitonic resonances mediated by non-adiabatic effects, which leads to a strong redistribution of scattering weight from the first to the second exciton, something that cannot be captured in an adiabatic theory. In the case of monolayer MoS2, we investigate (i) the difference in the coupling of the three main excitonic resonances to the in-plane and out-out-plane Raman modes and (ii) the effects of non-adiabaticity on the position of the main resonances.

Methods

Our computational approach to resonant Raman scattering is based on the general theoretical treatment we devised in Ref.  (?). In this work we are now interested in the special case of one-phonon Raman scattering. The rate for incident light of frequency $\omega_{\mathrm{L}}$ and polarization direction $\mu$ to scatter and be detected with frequency $\omega_{\mathrm{D}}$ and polarization $\nu$ after the creation of one phonon is given by:

[TABLE]

Here, the sum runs over all phonon modes $\lambda$ with zero momentum, as enforced by momentum conservation when the light is treated in the dipole approximation, in which the light momentum vanishes. Eq. (1) expresses the one-phonon part of the Raman spectrum as a sum of phonon spectral functions $\mathcal{A}_{\lambda}$, which determine the peak shapes, weighted by the scattering matrix element $\tilde{\mathcal{M}}^{\lambda}_{\mu\nu}(\omega_{\mathrm{L}}-\omega_{\mathrm{D}},\omega_{\mathrm{D}})$, which determines the peak intensities, or more precisely, their integrated areas. In this article, we focus on the scattering matrix element alone and treat the phonon spectral function as being a simple Lorentzian centered on the phonon frequency $\omega_{\lambda}$, with a width narrow enough to justify the approximation $\omega_{\mathrm{L}}=\omega_{\mathrm{D}}+\omega_{\lambda}$ in the following. The scattering matrix element contains information on the internal scattering dynamics and depends on the frequencies of both the incident and scattered light. This is crucial to correctly capture resonant scattering processes between the electronic excitations inside the material that mediate the scattering process. In Ref.  (?) we also gave a general expression for the scattering matrix element that mathematically relates it to the correlation of three observables: an electronic current generated by the incoming light, an electronic current that generates the outgoing light, and the electronic force on an atom that induces the lattice vibration and creates a phonon. More precisely, the scattering matrix element is given by the expression

[TABLE]

i.e., it is given by the double Fourier transform of a ”connected”, time-ordered correlation function of two zero-momentum electronic current densities $\hat{J}_{\mu,\nu}$ and an electronic force on the atoms, $\hat{F}_{\lambda}$ projected on the eigenmode of the phonon to be created. A precise definition of all the quantities that appear in it is given in Section S1 of the Supplementary Materials (SM).

The advantages of this approach to the description of Raman scattering are twofold. Firstly, it is a fully dynamical, i.e., it does not neglect any frequency dependence of the intensity as previous approaches do. This is important when the scattering dynamics involve scattering between electronic states that can be resonantly excited, not only with the light, but also with the phonon, which cannot be captured with existing, adiabatic approaches. Naively, these sort of non-adiabatic effects are argued to not be important in wide-gap semiconductors, as the size of the fundamental energy gap is much larger than typical phonon energies. However, beyond the fundamental gap, the validity of the adiabatic assumption should not be taken from granted. Indeed, as we show in the case of the hBN below, the phonon can resonantly couple two excited states, which has an a considerable effect on the Raman intensity. Secondly, our formulation permits a systematic approximation of the scattering matrix element, as the correlation function can be controllably approximated with methods from MBPT. This allows us to derive an expression for the scattering matrix element that includes both intermediate exciton states and their potentially resonant coupling by a phonon:

[TABLE]

This expression has the typical form of an expression that appears in third-order quantum mechanical perturbation theory and describes the excitation of an excitonic state $S$ with energy $E_{S}$ by the incoming light, its scattering to another excitonic state $S^{\prime}$ via emission of a phonon, and finally the recombination of exciton $S^{\prime}$ under emission of light to be detected. The exciton-light and exciton-phonon matrix elements in the numerators are defined as

[TABLE]

and

[TABLE]

respectively, where $\mathcal{C}$ ($\mathcal{V}$) denotes the set of all conduction (valence) band indices. $A^{S}_{\mathbf{k},c,v}$ denotes an excitonic eigenvector, corresponding to an envelope wave function in transition space, and $d^{\mu,(\mathrm{b})}_{\mathbf{k},a,b}$ and $g^{\lambda}_{\mathbf{k},a,b}$ denote the bare independent-particle dipole moment and screened electron-phonon coupling matrix elements, respectively. Finally, $\gamma$ is a broadening factor that we set to a constant value of $75$ meV, and the sums run over all excitonic states of positive energy. Details of how to obtain Eq. (3) from Eq. (2) are provided in Section S1 and Figs. S1 and S2 of the SM.

To make this approach computationally applicable, we implemented Eq. (3) in a Python code and obtain the necessary ingredients from first principles calculations. In general, we obtain the electronic band structure and wave functions on the level of DFT within the local density approximation in the parametrization due to Perdew and Zunger  (?). The DFT calculations were performed with the PWscf code of the QuantumESPRESSO suite  (?). Further material-specific computational details of the DFT calculations are provided in Section S2 of the SM.

The calculation of the excitonic eigenenergies and envelope wave functions was performed with the yambo code  (?). The two-particle interaction kernel was computed on the level of the Tamm-Dancoff approximation (TDA) as the sum of a statically screened, attractive Coulomb part plus a repulsive, bare exchange part. This approximation is routinely used in calculations of optical absorption spectra  (?, ?) and is sufficient to capture excitonic effects in hBN  (?), MoS2  (?), and many other materials. The static inverse dielectric screening function was computed on the level of the random phase approximation using plane-wave cutoffs of 60 Ry each for the electronic wave functions and the dielectric response function in reciprocal space. For the calculation of the two-particle interaction kernel, we truncated the attractive screened Coulomb part at a plane-wave energy cutoff of 10 Ry and the repulsive bare exchange part at a plane-wave cutoff of 20 Ry. We include a material-specific, finite number of bands in the interaction kernel and added a rigid energy shift to the independent-particle transition energies in order to account for the underestimation of the band gap in DFT. The material specific information is provided in Section S2 of the SM. Finally, we solve the Bethe-Salpeter equation iteratively for the first 2000 excitonic eigenenergies and envelope wave functions using the SLEPc library  (?).

For the screened electron-phonon coupling, we approximate the two-particle interaction kernel on the level of time-dependent density functional theory. As this approximation is equivalent to obtaining the screened electron-phonon matrix elements with density functional perturbation theory (see Secs. 4.3 and 4.5 of Ref.  (?)), we calculate the screened matrix elements with the PH code of the QuantumESPRESSO suite  (?). For the material specific phonon frequencies, we use the calculated values given in Section S2 of the SM.

Results and Discussion

As a first test case, we apply our approach to hexagonal boron nitride in its bulk form with $AA^{\prime}$ stacking order. This material features both strong excitonic effects  (?) and light atoms, which leads to relatively large phonon energies of up to $\approx$170 meV. The latter makes it an ideal candidate to observe the impact of non-adiabatic effects on the frequency dependence of the resonant Raman intensity.

In Fig. 1, we present the results of our calculation for the Raman-active, degenerate, in-plane optical ($E_{2\mathrm{g}}$) phonon mode, summed over outgoing photon and phonon polarizations and averaged over incoming photon polarizations. We show the computed Raman intensity as a function of excitation energy both for the full calculation with excitonic effects (blue lines) and for the independent-particle case (red lines), where we set the electron-hole interaction to zero while still retaining the TDA. In addition, we also compute the Raman intensity both in the full, non-adiabatic theory (full lines) and in the adiabatic limit (dashed lines). As previously shown  (?), the latter can be obtained from the former by neglecting the dependence of the scattering matrix element on the phonon frequency. For comparison, we also show the results of a calculation of the imaginary part of the dielectric function, which is closely related to the absorption coefficient. Additionally, we also included two experimental data points, measured in the pre-resonant regime, from Ref.  (?), which compare decently well to the results of our calculation.

In the IP-picture, the absorption spectrum of bulk hBN (red shade) sets on at the direct gap of around 6.50 eV located at the $M$-point in the band structure (see top panel of Fig. 2). When excitonic effects are included (blue shade), the absorption spectrum is dominated by a strong excitonic resonance near 5.58 eV, followed by a series of excitonic peaks with less oscillator strength (compare also Refs.  (?, ?, ?)).

By contrast, the Raman intensity shows a strikingly different behavior. On the IP-level, the Raman intensity rises only very slowly after the absorption onset. This behavior is typical for hexagonal materials and can be traced back to angular momentum conservation  (?, ?, ?). While the transitions at the direct gap at the $M$-point are relatively suppressed in the Raman process due to their small dipole moments (see top panel of Fig. 2), the energetically close transitions around the $H$-point cannot contribute for symmetry reasons. The latter transitions are located on a circularly symmetric part of the band structure, which implies full angular momentum conservation. This in turn suppresses the Raman scattering process for a degenerate phonon, as the latter carries an effective angular momentum of $\pm\hbar$ and cannot satisfy angular momentum conservation together with the incoming and outgoing photons. The same holds true for the transitions around the $K$-point, which also posses full rotation symmetry, with the exception of a small valley along the $\Gamma$-$K$ direction, which are the only transitions that can contribute to resonant low-energy Raman scattering on the IP-level.

When excitonic effects are included, we observe a strong suppression of the third and fourth excitonic peak. This suppression can be understood by analyzing their composition in $\mathbf{k}$-space (see bottom panel of Fig. 2). The envelope wave functions of these two excitons are strongly localized around the $K$- and $H$-points and are circularly symmetric with only weak trigonally warped contributions. As a result, they only weakly break the full rotation symmetry around the $K$-and $H$-points in the band structure and hence the Raman intensity for a degenerate phonon is again strongly suppressed at these frequencies due to angular momentum conservation. The first two excitons, by contrast, are more delocalized in $\mathbf{k}$-space and their wave functions show strong signs of trigonal warping. Thus the full rotation symmetry is broken down to the $120^{\circ}$-symmetry of the lattice, which then allows the Raman process on symmetry grounds.

Even more interesting though is the comparison to the results obtained within the adiabatic limit (dashed lines in Fig. 1). In this limit, the relative weight of the first two excitonic resonances follows that of the absorption spectrum. When non-adiabatic effects are included, however, we observe a strong redistribution of weight from the first to the second exciton. The reason for this is the non-negligible phonon energy of around 170 meV, which is close to the energy difference of the first two excitons and leads to strong, resonant inter-exciton scattering. In order to clearly demonstrate that the redistribution of oscillator strength from the second to the first bright exciton in the non-adiabatic theory indeed arises from inter-exciton scattering, we calculate the Raman intensity with only these contributions. For this, we define a Raman matrix element that only contains intra-exciton scattering processes via the restricted sum

[TABLE]

We then further define the Raman matrix element with only inter-exciton scattering processes as the difference

[TABLE]

of the full Raman matrix element of Eq. (3) and the intra-exciton-only one from Eq. (6). Note in this context that many of the exciton states are degenerate. For example, the first excitonic resonance in the absorption spectrum is due to a bright exciton doublet of $E_{2\mathrm{u}}$-symmetry, which is slightly higher in energy than the lowest (dark) exciton doublet of $E_{2\mathrm{g}}$-symmetry  (?). The double sum over excitons in Eq. (6) is therefore understood to run over all exciton states, taking into account all off-diagonal exciton-phonon matrix elements between degenerate states, but ignoring the coupling of one exciton multiplet to another.

With these two restricted matrix elements, we can compute the Raman intensity again, this time, however, only with intra- or inter-exciton scattering. The resulting intensities are shown in Fig. 3. Firstly, it is important to note that the total intensity is not equal to the sum of intensities with only one of the two scattering channels, as the two contributions are first added to the total matrix element and then squared. In the case of the adiabatic limit (bottom panel), intra- (lighter color) and inter-exciton scattering (darker color) contribute to the total Raman intensity with about equal weight and they add up constructively. By contrast, the full non-adiabatic result (top panel) is strongly dominated by the contribution from *inter-*exciton scattering. Moreover, while at the first excitonic resonance, inter- and intra-exciton scattering reinforce each other, at the second excitonic resonance, the two contributions destructively interfere. Compared to the adiabatic case, we thus observe that the non-adiabatic contribution from the inter-exciton scattering terms changes sign, which is a typical behavior around a resonance and thus further demonstrates the importance of considering resonant scattering due to the finite phonon frequency. This illustrates that non-adiabatic effects can play a significant role even in wide-band gap materials and emphasizes the importance of considering both non-adiabatic and excitonic effects at the same time.

However, there are also cases in which the resonant, non-adiabatic coupling can be suppressed by other effects. One example is the case of monolayer molybdenum disulfide which is also known to display strong excitonic effects  (?), yet in which spin-orbit coupling also plays a significant role. MoS2 features two Raman-active modes, a degenerate in-plane mode of symmetry $E^{\prime}$ and an out-of-plane one of symmetry $A^{\prime}_{1}$. Furthermore, MoS2 possesses an electronic band structure that in its optically active region is strongly influenced by spin-orbit coupling  (?). In Fig. 4, we show the results of our calculation for both Raman-active modes, with computational details of the underlying ab initio calculations provided in Section S2 of the SM. The experimental data shown in the figure was taken from Ref.  (?).

The optical absorption spectrum of monolayer MoS2 (blue shade in Fig. 4) is dominated by three strong excitonic resonances, at around 1.97 eV, 2.11 eV, and in the region 2.6–2.8 eV. The first two of these, which each are doubly degenerate, are commonly referred to as the “A”- and “B”-exciton  (?, ?), while the latter resonance actually consists of several excitons and which are collectively referred to as the “C”-exciton. The $\mathbf{k}$-point-resolved envelope wave functions of the “A”- and “B”-excitons are strongly localized around the $K$-point (see bottom panel of Fig. 5) and the difference in energy corresponds to the spin-orbit coupling-induced splitting of the valence band. By contrast, the dominant “C”-exciton is strongly delocalized in $\mathbf{k}$-space. Note that in our calculations, the onset of the “C”-exciton was obtained at too high an energy. In order to facilitate a comparison of the relative heights of the Raman peaks, we compressed the excitonic spectrum in excess of 2.25 eV to match the experimental peak positions in absorption spectra.

The Raman spectrum of the $A^{\prime}_{1}$-mode (blue line) closely follows the trend of the absorption spectrum. Due to the double-resonant structure of the scattering matrix element, the excitonic peaks become more pronounced and well-separated, for example, in the case of the “A’- and “B”-excitons. The degenerate $E^{\prime}$-mode, on the other hand, does not show any resonances at the energies of the “A’- and “B”-excitons. This can again be traced back to their envelope wave functions being strongly localized around the $K$-point, where the transition band structure exhibits full rotation symmetry (see top panel of Fig. 5), which in turn suppresses Raman scattering with the $E^{\prime}$-mode due to angular momentum conservation. This trend is also seen in the experimental data, in which the $A^{\prime}_{1}$-mode features much more prominently in the spectrum. Note that the peak positions in our calculations appear to be shifted by around 80 meV. One reason for this could be the missing inclusion of higher-order electron-phonon coupling effects, such as the zero-point vibrational renormalization of the electronic band structure, which is known to reduce the fundamental band gap by around 75 meV  (?). In the region around the “C”-exciton, however, the excitons couple much more strongly to the $E^{\prime}$-mode than to the $A^{\prime}_{1}$-mode. This behavior can again be understood from an analysis in terms of inter- and intra-exciton scattering (see Fig. 6).

In Fig. 6, we show the Raman intensity with only one of the two contributions for both the $E^{\prime}$-mode (top panel) and the $A^{\prime}_{1}$-mode (bottom mode). We only show the results for the non-adiabatic theory, as the adiabatic limit yields almost the same results.

Firstly, we note that in case of the $E^{\prime}$-mode, the Raman intensity is strongly suppressed around the “A”- and “B”-excitons, irrespective of whether only intra- or inter-exciton scattering is considered alone. This further underlines that the suppression is indeed the result of a physical symmetry, i.e., in this case approximate full angular momentum conservation. Secondly, the Raman intensity for the $E^{\prime}$-mode in the energy regime close to the “C”-exciton is dominated inter-exciton scattering and intra-exciton scattering plays almost no role. In the case of the $A^{\prime}_{1}$-mode, the picture is entirely different, however. The low-energy regime around the “A”- and “B”-excitons is governed by intra-exciton scattering due to the isolated nature of the two excitons. At higher excitation energies, inter-exciton scattering gains in weight and becomes almost as important as intra-exciton processes. Most importantly though is the fact that the two scattering mechanism contribute to the total Raman matrix element with opposite signs. As such, the lower intensity of the $A^{\prime}_{1}$-mode compared to that of the $E^{\prime}$-mode is the effect of destructive quantum interference.

This results in the $E^{\prime}$-mode featuring a much larger intensity than the $A^{\prime}_{1}$-mode in the exciton continuum. By contrast, the isolated nature of the “A”- and “B”-excitons means that inter-exciton scattering plays no role and only scattering within either the “A”- or “B”-exciton doublet is important. This analysis also corrects the controversial statement of the recent Ref.  (?) that only inter-exciton scattering should be considered in the resonant regime, i.e., the “A”- and “B”-excitons should be silent. This is not only in contrast to the analysis presented here, but is also in stark contrast to experiments of other monolayer transition metal dichalcogenites  (?).

Finally, we comment on the effects of non-adiabaticity in the case of MoS2. Here, the phonon energies of the Raman active modes are on the order of 50 meV  (?). This is not enough to overcome the spin-orbit coupling-induced splitting of the “A”- and “B”-excitons of around 110 meV. The contributions from the “A”- and “B”-excitonic resonances to the scattering matrix element are thus prevented from interfering with one another.

To verify this hypothesis, we artificially triple the phonon frequency to overcome the gap between the “A”-to-“B” excitons and computed the Raman intensity again (see purple line in inset in bottom panel of Fig. 6). Now the intensity decreases significantly as a result of destructive interference between the two excitonic resonances, which are now resonantly coupled. We can thus conclude that in the real case, the spin-orbit splitting suppresses the resonant coupling of the two excitons with a phonon, leading to a bigger Raman intensity. As such, the inclusion of non-adiabatic effects only results in a slight blueshift of all resonances compared to the adiabatic case (dashed lines in Fig. 4).

In conclusion, we have presented a method for the computation of resonant one-phonon Raman intensities from first principles that takes into account both excitonic and non-adiabatic effects. We have applied our method to both bulk hexagonal boron nitride and single-layer molybdenum disulfide. In the case of the former, we have explained the absence of particular higher excitonic resonances from the Raman intensity spectrum with angular momentum conservation. More importantly, we have proven the significance of non-adiabaticity for the Raman process, as it leads to strong quantum interference between the first two excitonic resonances. For MoS2, we have shown that spin-orbit splitting protects the lowest two excitonic resonances from non-adiabatic effects and that the in-plane $E^{\prime}$-mode is silent in this regime due to angular momentum conservation. Furthermore, we have explained the different behavior of the two Raman active modes at higher excitation energies by destructive quantum interference of intra- and inter-exciton scattering processes involving the out-of-plane $A^{\prime}_{1}$-mode.

Acknowledgments

The authors acknowledge the use of the QuantumESPRESSO suite  (?) for the DFT calculations within the local density approximation  (?), the yambo code  (?) and the SLEPc  (?) library for the MBPT calculations, and usage of the HPC facilities of the University of Luxembourg  (?). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data available from authors upon request. The authors would like to thank A. Marini for inspiring discussions. S.R. and L.W. acknowledge financial support by the National Research Fund (FNR) Luxembourg (projects RAMGRASEA and INTER/ANR/13/20/NANOTMD). The authors declare that they have no competing interests. Author contributions: S.R. and L.W. devised the project idea. S.R. developed the theory, implemented it, performed the calculations, and prepared the main part of the manuscript. Both authors discussed the results, the ideas for their analysis, and edited the manuscript.

Supplementary materials

The supplementary materials include the following items:

Section S1 – Theoretical details

Figure S1 – Diagrammatic analysis of the scattering matrix element

Figure S2 – Ladder-like approximation to the scattering matrix element

Section S2 – Computational details for the hBN and MoS2 calculations

See pages - of supplementary