Research ArticlePHYSICS

Nonadiabatic exciton-phonon coupling in Raman spectroscopy of layered materials

See allHide authors and affiliations

Science Advances  07 Aug 2020:
Vol. 6, no. 32, eabb5915
DOI: 10.1126/sciadv.abb5915


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 MoS2, we observe that quantum interference effects are suppressed by the spin-orbit splitting of the excitons.


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 (713), 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, 1418), 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 (MoS2), 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 MoS2, 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.


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Ṗ1ph.=const.×ωLωDλM˜μνλ(ωLωD,ωD)2Aλ(ωLωD)(1)

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 Mμνλ(ωLωD,ωD), which determines the peak intensities or, more precisely, their integrated areas. Here, we focus on the scattering matrix element alone and treat the phonon spectral function as being a simple Lorentzian centered on the phonon frequency ωλ, 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 expressionM˜μνλ(ω,ω)1ħ2+dt eiωt+dt eiωt(i)0T[Fˆλ(t)Jˆν(t)Jˆμ(0)]0connect.(2)i.e., it is given by the double Fourier transform of a “connected,” time-ordered correlation function of two zero-momentum electronic current densities Ĵμ,ν and an electronic force on the atoms, F̂λ, 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.

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 phononM˜μνλ(ωλ,ωLωλ)=S,SʹdSμ(gS,Sʹλ)*(dSʹν)*(ħωLES+iγ)(ħωLħωλESʹ+iγ)+S,Sʹ(dSμ)*gS,SʹλdSʹν(ħωL+ESiγ)(ħωLħωλ+ESʹiγ)(3)

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 ES 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 asdSμk,bVaCdk,a,bμ,(b)(Ak,a,bS)*(4)andgS,Sλk[a,cCbV(Ak,a,bS)*gk,a,cλAk,c,bSaCb,cV(Ak,a,bS)*gk,c,bλAk,a,cS](5)respectively, where C (V) denotes the set of all conduction (valence) band indices. Ak,c,vS denotes an excitonic eigenvector, corresponding to an envelope wave function in transition space, and dk,a,bμ,(b) and gk,a,bλ denote the bare independent-particle (IP) dipole moment and screened electron-phonon coupling matrix elements, respectively. Lastly, γ 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.

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), MoS2 (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.


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 (E2g) 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.

Fig. 1 Calculated Raman intensities for hBN as a function of excitation frequency.

Blue lines, Raman intensity on the BSE level calculated with the full nonadiabatic theory (full lines) and in the adiabatic (Ad.) limit (dashed lines). Red lines, same as blue lines, but without excitonic effects, i.e., on the IP level. Shaded areas, imaginary part of the dielectric function on the BSE (blue) and IP (red) levels. Squares, experimental (Exp.) data from (22), normalized such that the last point matches the theoretical curve. a.u., arbitrary units.

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 (3032)].

Fig. 2 hBN: Electronic transition band structure and exciton envelope wave functions in reciprocal space.

Top: Electronic transition band structure obtained from Kohn-Sham DFT within the LDA. A rigid shift of 2 eV was added to all conduction band energies. The opacity of the lines corresponds to the optical activity of the transition. Bottom: Exciton envelope wave functions in reciprocal space for the first four bright excitons. Each panel shows the band-summed square of the envelope wave function Ak,c,vS2 as a function of the wave vector k integrated over the kz coordinate. coord., coordinate.

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 ±ħ 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 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 sumM˜μνλintra(ωλ,ωLωλ)=S,SES=ES[dSμ(gS,Sλ)*(dSν)*(ħωLES+iγ)(ħωLħωλESʹ+iγ)+(dSμ)*gS,SʹλdSʹν(ħωL+ESiγ)(ħωLħωλ+ESʹiγ)](6)

We then further define the Raman matrix element with only inter-exciton scattering processes as the differenceM˜μνλinter(ωλ,ωLωλ)M˜μνλ(ωλ,ωLωλ)M˜μνλintra(ωλ,ωLωλ)(7)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 E2u symmetry, which is slightly higher in energy than the lowest (dark) exciton doublet of E2g 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.

Fig. 3 hBN: Raman intensity with only inter- or intra-exciton scattering.

Raman intensity with all contributions (thick lines) and with only inter- or intra-exciton scattering (darker and lighter thin lines, respectively) taken into account. Top: Full, nonadiabatic case. Bottom: Adiabatic limit. The shaded area represents the absorption spectrum. The Raman intensity is given in the same units in both panels and as in Fig. 1.

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. MoS2 features two Raman-active modes, a degenerate in-plane mode of symmetry E′ and an out-of-plane one of symmetry A1. Furthermore, MoS2 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).

Fig. 4 Calculated Raman intensities for monolayer molybdenum disulfide as a function of excitation frequency.

Green and blue lines, Raman intensity of the E′-mode (green) and A1-mode (blue) on the BSE level calculated with the full nonadiabatic theory (full lines) and in the adiabatic limit (dashed lines). Blue shaded area, imaginary part of the dielectric function on the BSE level. Green squares and blue triangles, experimental data from (5), normalized such that the A1-mode intensity at the second (“B”)-exciton matches the theoretical curve. Inset: Zoom-in into the region around the first two excitons, highlighting an 80-meV shift between theory and experiment.

The optical absorption spectrum of monolayer MoS2 (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.

Fig. 5 MoS2: Electronic transition band structure and exciton envelope wave functions in reciprocal space.

Top: Electronic transition band structure obtained from Kohn-Sham DFT within the LDA. A rigid shift of 0.925 eV was added to all conduction band energies. The opacity of the lines corresponds to the optical activity of the transition. Bottom: Exciton envelope wave functions in reciprocal space for the first three bright excitons (the A-, B-, and C-excitons). Each panel shows the band-summed square of the envelope wave function Ak,c,vS2 as a function of the wave vector k. In the case of the C-exciton, we show the sum of the squared wave functions of five energetically close bright excitons.

The Raman spectrum of the A1-mode (blue line) closely follows the trend of the absorption spectrum. Because of 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′-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 A1-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 bandgap by around 75 meV (21). In the region around the C-exciton, however, the excitons couple much more strongly to the E′-mode than to the A1-mode. This behavior can again be understood from an analysis in terms of inter- and intra-exciton scattering (see Fig. 6).

Fig. 6 MoS2: Raman intensity with only inter- or intra-exciton scattering.

Raman intensity with all contributions (thick lines) and with only inter- or intra-exciton scattering (darker and lighter thin lines, respectively) taken into account. Top: E′-mode. Bottom: A1-mode. The shaded area represents the absorption spectrum. The Raman intensity is given in the same units in both panels and as in Fig. 4. Insets: Zoom-ins into the areas around the three dominant excitonic resonances. In the bottom inset, the purple line depicts the Raman intensity of the A1-mode with an artificially tripled phonon frequency (ph. freq.).

In Fig. 6, we show the Raman intensity with only one of the two contributions for both the E′-mode (top) and the A1-mode (bottom). We only show the results for the nonadiabatic theory, as the adiabatic limit yields almost the same results.

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 A1-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 mechanisms contribute to the total Raman matrix element with opposite signs. Hence, the lower intensity of the A1-mode compared to that of the E′-mode is the effect of destructive quantum interference.

This results in the E′-mode featuring a much larger intensity than the A1-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 study (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 MoS2. 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 MoS2, 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 A1-mode.


Supplementary material for this article is available at

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.


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.

Stay Connected to Science Advances

Navigate This Article