## Abstract

We present an ab initio computational approach for the calculation of resonant Raman intensities, including both excitonic and nonadiabatic 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.

## INTRODUCTION

Raman spectroscopy is a versatile tool for the characterization of many materials and, in particular, two-dimensional materials (*1*). It allows the study of both vibrational and electronic properties. Among many other things, it serves to probe many-body (*2*) and nonadiabatic effects (*3*, *4*) as well as excitonic ones (*5*, *6*). 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 (*7*–*13*), which is able to capture nonadiabatic 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 (*6*, *14*–*18*), which captures excitonic effects but is entirely adiabatic. Hence, an approach that is able to capture both effects is highly desirable not only from a theoretical and conceptual point of view but also in practice. This is especially true for systems with sizable phonon frequencies and strong excitonic effects and also for materials whose properties are a priori unknown.

In a recent work (*19*), 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 nonadiabatic 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 quasiparticle band structure and solving the Bethe-Salpeter equation (BSE) to include excitonic effects. We then apply it to two materials, bulk hexagonal boron nitride (hBN) and monolayer molybdenum disulfide (MoS_{2}), which are known to feature strong excitonic effects in both absorption (*20*, *21*) and Raman spectra (*5*, *22*) 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 nonadiabatic 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 MoS_{2}, we investigate (i) the difference in the coupling of the three main excitonic resonances to the in-plane and out-of-plane Raman modes and (ii) the effects of nonadiabaticity on the position of the main resonances.

## METHODS

Our computational approach to resonant Raman scattering is based on the general theoretical treatment that we devised in (*19*). In this work, we are now interested in the special case of one-phonon Raman scattering. The rate for incident light of frequency ω_{L} and polarization direction μ to scatter and be detected with frequency ω_{D} and polarization ν after the creation of one phonon is given by

Here, the sum runs over all phonon modes λ with zero momentum, as enforced by momentum conservation when the light is treated in the dipole approximation, in which the light momentum vanishes. Equation 1 expresses the one-phonon part of the Raman spectrum as a sum of phonon spectral functions A_{λ}, which determine the peak shapes, weighted by the scattering matrix element _{λ}, with a width narrow enough to justify the approximation ω_{L} = ω_{D} + ω_{λ} 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 (*19*), 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

The advantages of this approach to the description of Raman scattering are twofold. First, it is 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 nonadiabatic 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 for granted. As we show in the case of hBN below, the phonon can resonantly couple two excited states, which has a considerable effect on the Raman intensity. Second, 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

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*′ via emission of a phonon, and lastly the recombination of exciton

*S*′ under emission of the light to be detected. The exciton-light and exciton-phonon matrix elements in the numerators are defined as

To make this approach computationally applicable, we implemented Eq. 3 in a `Python` code and obtained 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 (LDA) in the parameterization due to Perdew and Zunger (*23*). The DFT calculations were performed with the `PWscf` code of the `Quantum ESPRESSO` suite (*24*). Further material-specific computational details of the DFT calculations are provided in section S2.

The calculation of the excitonic eigenenergies and envelope wave functions was performed with the `yambo` code (*25*). 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 (*26*, *27*) and is sufficient to capture excitonic effects in hBN (*20*), MoS_{2} (*21*), 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 Rydberg (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 to account for the underestimation of the bandgap in DFT. The material-specific information is provided in section S2. Last, we iteratively solve the Bethe-Salpeter equation for the first 2000 excitonic eigenenergies and envelope wave functions using the `SLEPc` library (*28*).

For the screened electron-phonon coupling, we approximate the two-particle interaction kernel on the level of time-dependent DFT. As this approximation is equivalent to obtaining the screened electron-phonon matrix elements with density functional perturbation theory [see sections 4.3 and 4.5 of (*29*)], we calculate the screened matrix elements with the `PH` code of the `Quantum ESPRESSO` suite (*24*). For the material-specific phonon frequencies, we use the calculated values given in section S2.

## RESULTS AND DISCUSSION

As a first test case, we apply our approach to hBN in its bulk form with *AA*′ stacking order. This material features both strong excitonic effects (*20*) and light atoms, which leads to relatively large phonon energies of up to ≈170 meV. The latter makes it an ideal candidate to observe the impact of nonadiabatic 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*_{2g}) 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, nonadiabatic theory (full lines) and in the adiabatic limit (dashed lines). As previously shown (*19*), 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. In addition, we also included two experimental data points, measured in the preresonant regime, from (*22*), 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 (*30*–*32*)].

By contrast, the Raman intensity shows a notably 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 (*6*, *13*, *33*). While the transitions at the direct gap at the *M* point are relatively suppressed in the Raman process because of 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 *K* point, which also have full rotation symmetry, with the exception of a small valley along the Γ − *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 peaks. This suppression can be understood by analyzing their composition in *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 because of angular momentum conservation. The first two excitons, by contrast, are more delocalized in *k*-space, and their wave functions show strong signs of trigonal warping. Thus, the full rotation symmetry is broken down to the 120^{∘} 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 nonadiabatic 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. To clearly demonstrate that the redistribution of oscillator strength from the second to the first bright exciton in the nonadiabatic 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

We then further define the Raman matrix element with only inter-exciton scattering processes as the difference*E*_{2u} symmetry, which is slightly higher in energy than the lowest (dark) exciton doublet of *E*_{2g} symmetry (*32*). 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. First, 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 (Fig. 3, bottom), intra-exciton scattering (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 nonadiabatic result (Fig. 3, top) 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, and at the second excitonic resonance, the two contributions destructively interfere. Compared to the adiabatic case, we thus observe that the nonadiabatic 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 nonadiabatic effects can play a crucial role even in wide-bandgap materials and emphasizes the importance of considering both nonadiabatic and excitonic effects at the same time.

However, there are also cases in which the resonant, nonadiabatic 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 (*31*), yet in which spin-orbit coupling also plays a crucial role. MoS_{2} features two Raman-active modes, a degenerate in-plane mode of symmetry *E*′ and an out-of-plane one of symmetry _{2} has an electronic band structure that, in its optically active region, is strongly influenced by spin-orbit coupling (*21*). 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. The experimental data shown in the figure were taken from (*5*).

The optical absorption spectrum of monolayer MoS_{2} (blue shade in Fig. 4) is dominated by three strong excitonic resonances, at around 1.97 and 2.11 eV and in the region between 2.6 and 2.8 eV. The first two of these, which each are doubly degenerate, are commonly referred to as the “A”- and “B”-exciton (*34*, *35*), while the latter resonance actually consists of several excitons, which are collectively referred to as the “C”-exciton. The *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 *k*-space. Note that in our calculations, the onset of the C-exciton was obtained at too high an energy. 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 *E*′-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*′-mode due to angular momentum conservation. This trend is also seen in the experimental data, in which the *21*). In the region around the C-exciton, however, the excitons couple much more strongly to the *E*′-mode than to the

In Fig. 6, we show the Raman intensity with only one of the two contributions for both the *E*′-mode (top) and the

First, we note that in the case of the *E*′-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. Second, the Raman intensity for the *E*′-mode in the energy regime close to the C-exciton is dominated by inter-exciton scattering, and intra-exciton scattering plays almost no role. In the case of the *E*′-mode is the effect of destructive quantum interference.

This results in the *E*′-mode featuring a much larger intensity than the *18*) 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 also in stark contrast to experiments of other monolayer transition metal dichalcogenides (*6*).

Lastly, we comment on the effects of nonadiabaticity in the case of MoS_{2}. Here, the phonon energies of the Raman-active modes are on the order of 50 meV (*36*). 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 and B excitons and computed the Raman intensity again (see purple line in inset in the 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. Hence, the inclusion of nonadiabatic 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 nonadiabatic effects. We have applied our method to both hBN 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. We have proven the crucial role of nonadiabaticity for the Raman process, as it leads to strong quantum interference between the first two excitonic resonances. For MoS_{2}, we have shown that spin-orbit splitting protects the lowest two excitonic resonances from nonadiabatic effects and that the in-plane *E*′-mode is silent in this regime because of 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

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/32/eabb5915/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We acknowledge the use of the

`Quantum ESPRESSO`suite (

*24*) for the DFT calculations within the local density approximation (

*23*), the

`yambo`code (

*25*), and the

`SLEPc`library (

*28*) for the MBPT calculations and usage of the HPC facilities of the University of Luxembourg (

*37*). We would like to thank A. Marini for inspiring discussions.

**Funding:**S.R. and L.W. acknowledge financial support from the National Research Fund (FNR) Luxembourg (projects RAMGRASEA and INTER/ANR/13/20/NANOTMD).

**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 and the ideas for their analysis and edited the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).