Ab Initio Approach to Second-order Resonant Raman Scattering Including Exciton-Phonon Interaction
Yannick Gillet, Stefan Kontur, Matteo Giantomassi, Claudia Draxl, Xavier Gonze
aa r X i v : . [ c ond - m a t . m t r l - s c i ] A ug Ab Initio
Approach to Second-order ResonantRaman Scattering Including Exciton-PhononInteraction
Yannick Gillet , Stefan Kontur , Matteo Giantomassi , Claudia Draxl , and XavierGonze Universit ´e catholique de Louvain, Institute of Condensed Matter and Nanosciences, Nanoscopic Physics,Chemin des ´etoiles 8, bte L7.03.01, 1348 Louvain-la-Neuve, Belgium Physics Department and IRIS Adlershof, Humboldt-Universit ¨at zu Berlin, Zum Großen Windkanal 6, D-12489Berlin, Germany European Theoretical Spectroscopy Facility (ETSF) * [email protected] ABSTRACT
Raman spectra obtained by the inelastic scattering of light by crystalline solids contain contributions from first-order vibrationalprocesses (e.g. the emission or absorption of one phonon, a quantum of vibration) as well as higher-order processes with atleast two phonons being involved. At second order, coupling with the entire phonon spectrum induces a response that maystrongly depend on the excitation energy, and reflects complex processes more difficult to interpret. In particular, excitons(i.e. bound electron-hole pairs) may enhance the absorption and emission of light, and couple strongly with phonons inresonance conditions. We design and implement a first-principles methodology to compute second-order Raman scattering,incorporating dielectric responses and phonon eigenstates obtained from density-functional theory and many-body theory.We demonstrate our approach for the case of silicon, relating frequency-dependent relative Raman intensities, that are inexcellent agreement with experiment, to different vibrations and regions of the Brillouin zone. We show that exciton-phononcoupling, computed from first principles, indeed strongly affect the spectrum in resonance conditions. The ability to analyzesecond-order Raman spectra thus provides direct insight into this interaction.
Introduction
In Raman spectroscopy, the change of frequency and intensity of light, that is inelastically scattered by the material, allowsone to collect a rich set of information, even more if the effect is monitored as a function of the frequency of the incominglight . In most systems, the energy loss is due to the creation of one or more phonons, and different orders of phonon-photon scattering processes contribute to the scattering cross section. The first-order contribution, in which only one phonon isemitted or absorbed, yields information on vibrations with very low crystalline momentum, due to the negligible momentumof light. Extracting information on phonon frequencies from the first-order spectrum is straightforward, as clear selectionrules govern the appearance or extinction of usually well-isolated peaks. In contrast, second-order Raman processes exhibitfeatures coming from different crystalline momenta, potentially from the entire Brillouin zone (BZ), with the only constraintof negligible momentum of the phonon pairs involved in the two-phonon process. They are more complex and, hence, moredifficult to interpret. When the frequency of the incoming light comes close to the specific frequency needed to drive thetransfer of an electron from an occupied state to an unoccupied state, which is termed a resonance, the absolute Ramanintensities can change by several orders of magnitude. This resonance phenomenon appears independently of the numberof emitted phonons. Thus a particular frequency can be selected to increase the Raman signal . The possibility to createbound electron-hole pairs, called excitons, crucially impact optical absorption and related photocurrent in many materials, amechanism at work in commercial photovoltaic devices . Understanding and characterizing their dynamics and interactionsis an experimental challenge. In this context, the ability to interpret the widely available Raman data to characterize exciton-phonon interaction is crucial for the understanding of light-matter interaction in general, and in particular in technologicallyrelevant materials.Experimentally, the second-order part of the Raman spectrum is often used to characterize low-dimensional or layered ma-terials, such as e.g. graphene or graphite , and more recently also transition-metal dichalcogenides . The few existingfirst-principles studies for these systems relied on the independent-electron approximation (IPA) , in which the formationof excitons is neglected. As these low-dimensional materials exhibit very well-defined and specific spectral features , their1econd-order Raman spectra can still be interpreted quite easily, at least concerning the location of the peaks. However, forphotovoltaic applications, bulk materials form the vast majority of relevant materials. Their Raman spectra contains a largernumber of features, and is thus difficult to decipher. Resonant first-order and second-order Raman scattering processes mayhave e.g. different temperature behavior , whose most natural interpretation invoke characteristics of the direct and indirectgaps. For classical semiconductors such as silicon (and Ge, III-V etc), the analysis of second-order Raman spectra was facil-itated by the use of simple Hamiltonians, based on the effective mass approximation or on semi-empirical pseudopotentials,with model phonon band structure, model electron-phonon interaction and, in some cases, the inclusion of Wannier-Mott exci-tons . Early first-principles calculations were done for vanishing laser frequencies (no resonance effects) without excitoniceffects . Our work allows one to compute from first-principles second-order Raman spectra without any of these limitations.A prototypical example of resonant second-order Raman spectra is presented in Fig. 1. It shows the experimentalfrequency-dependent Raman spectra of silicon between 900 and 1050 cm − , as measured by Renucci and coworkers . Thespectra are characterized by significant variations in intensity depending on the frequency of the incoming light, even for fre-quencies corresponding to energies smaller than the direct gap of 3.2 eV. Such resonance effects are allowed in second-orderRaman processes, due to indirect transitions between electronic states that are caused by phonons of finite wavevector. Hence,it is clear that multi-phonon processes beyond first order are crucial. Below, we will show that all the features can be wellinterpreted by our first-principle theoretical approach.
900 950 1000 1050 1100 900 950 1000 1050 1100Raman shift (cm − )3.41 eV2.81 eV2.6 eV2.41 eV2.18 eV1.92 eV1.65 eV1.17 eV BSEIPAExperimentRaman shift (cm − ) I n t e n s i t y ( a r b . un i t s ) Figure 1. (Color online). Experimental (left, taken from Ref. 26) and theoretical (right) second-order Raman spectra ofsilicon, between 900 and 1050 cm − , for different excitation energies between 1.17 eV and 3.41 eV, normalized to theirmaximal values. For clarity, the curves are shifted with respect to each other (lowest excitation energy at the bottom).Theoretical results have been obtained in the independent-particle approximation (IPA) and in the Bethe-Salpeter equation(BSE) formalism.A general methodology for describing multi-phonon scattering processes up to arbitrary order has been proposed by Knolland Draxl , but its application had been long restricted to the ab initio calculations of first-order resonant Raman scattering, inthe IPA, without excitonic effects, e.g. in various superconducting materials and ladder compounds . Recently, thanksto the treatment of excitons using the Bethe-Salpeter equation (BSE), excitonic effects were included from first principles insuch first-order resonant Raman scattering . Excitonic effects were clearly observed in the absolute intensity enhancement,and proved crucial to obtain good agreement with experimental data. Other first-principles studies of Raman intensities havebeen based on density-functional perturbation theory . However, this last approach has been developed only in the limitof vanishing light frequency, so without addressing any resonance effect and without excitonic effects .In this work, we address the even more challenging computation of second-order resonant Raman scattering, includingexciton-phonon coupling. Even by today standards, the computational load is quite formidable. We compute the ingredientsusing state-of-the-art methodology (density-functional perturbation theory and Bether-Salpeter equation), while also showingthe adequacy of the simplifying “overtone” approximation. Applying it to the case of Si, we not only reproduce the abovedescribed experimental data, but give insight into the origin of the spectral features, and explore the role of excitonic effectson the Raman intensities. .0 1.5 2.0 2.5 3.010 BSEIPA1.0 2.0 3.00246 I
BSE /I IPA
Excitation energy ~ ω L (eV) I n t e g r a t e d i n t e n s i t y ( a r b . un i t s ) Figure 2. (Color online). Integrated intensities of BSE and IPA spectra over Raman shifts ranging from 900 cm − to 1050cm − , as a function of the excitation energy ¯ h ω L . The ( ω L − ω R ) prefactor, see Eq. 1 has not been included, and thespectrum has been renormalized with respect to the value at ¯ h ω L = . Results
We first describe the normalized intensities for Raman shifts in the 900-1100 cm − range. Such results are presented in theright panel of Fig. 1. On both levels of theory, IPA and BSE, the evolution of the intensity of the Raman shift with increasinglight frequency (or photon excitation energy) follows very closely the experimental data. While most of the spectral weightis localized around 950 cm − for the lowest excitation energy of 1.17 eV, the shoulder below 1000 cm − is getting morepronounced with increasing excitation energy and starts to dominate above 2.4 eV. A third peak around 1030 cm − that ishardly visible in the bottom curves, is particularly pronounced only for the highest excitation energy of 3.41 eV. In excellentagreement with experiment, its intensity is nearly the same as the peak at 950 cm − , and somewhat lower than the one at 980cm − .Second, we examine the integrated intensities. Fig. 2 shows a comparison between the IPA and the BSE absolute intensity,divided by ( ω L − ω R ) (see Methods section for more detail), integrated over the range between 900 and 1100 cm − , as afunction of the excitation energy. Even if in this range of frequencies the relative intensities computed in the two formalismsare very similar as shown in Fig. 1, their absolute intensities differ by a factor 5 for excitation energies around 3 eV, asexhibited in the inset of the Fig. 2. Clearly, excitonic effects are important, which will be analyzed in the Discussion section.The Raman spectrum in a broader range of Raman shifts is also obtained. In addition to the enhancement of the absoluteintensity, similar to the one shown in Fig. 2, Fig. 3 shows that the use of the BSE, that is, the inclusion of exciton-phononcoupling, results in a pronounced intensity increase in the part of the Raman spectrum between 600 and 900 cm − , especiallyfor excitation energies between 2 and 3 eV. This additional excitonic effect will also be discussed later.In order to compute these results, we have made the approximation that the second-order contribution to Raman spectrum isdominated by the emission of two phonons that are connected by time-reversal symmetry: same frequency, and time-reversedcollective eigendisplacement and momentum. From a group-theoretical point of view, derivatives with respect to such phononpairs will always be active and contribute to the fully symmetric component of the spectrum . However, group-theoreticalarguments do not insure the vanishing of other contributions for arbitrary phonon wavevectors. While this approximation isjustified a posteriori by the agreement show in Fig. 1, it has been nevertheless possible to theoretically examine its validity,in the specific case of vanishing excitation energy. We tested this approximation by computing the second-order derivativesof the dielectric constant obtained within Density-Functional Perturbation Theory (DFPT), where the electric field is usedas a perturbation . In this approach, the excitation energy vanishes, but unlike in the IPA, the local-field effects are takeninto account, which makes its level of sophistication intermediate between IPA and BSE. Fig. 4 shows the dominance of thederivatives by the overtone pair of modes (by around one order of magnitude with respect to other contributions in the worstcase), which fully justifies our approach. This ’overtone’ approximation was a crucial step in making the BSE computationtractable, as discussed in the Methods section, and to be able to treat excitonic effects. The confirmation of the validity ofthe ’overtone’ approximation for silicon, discussed in prior works, also for other simple semiconductors isthus another significant result of our work. However, there are several cases in which features can be assigned to combination .6 eV2.41 eV2.18 eV1.92 eV1.65 eV1.17 eV200 1100500 800 I n t e n s i t y ( a r b . un i t s ) Raman shift (cm − )BSEIPA Figure 3. (Color online). Second-order Raman spectra calculated within the IPA (top) and the BSE (bottom) for severalexcitation energies ¯ h ω L . Each spectrum is normalized with respect to its highest value, as in Fig. 1.
200 300 400 500 600 700 800 900 1000 1100 I n t e n s i t y ( a r b . un i t s ) Raman shift (cm − )OvertoneFull Figure 4. (Color online). Second-order Raman spectrum obtained by differentiating the dielectric constant (at ¯ h ω L = odes, like in Refs. 16, 47, 48. In such case, the identification of relevant combinations might be done in the static limit usingDFPT, followed by the computation of the contributions from specific mode combinations in the relevant energy spectrum (asis done already in our work, for degenerate phonon modes). Discussion
Our first-principles theory can provide insight into the origin of the features seen in Fig. 1, and then can shed light on theunderlying physics, especially the role of excitonic effects. For this purpose, we first present the expression that yields thetheoretical Raman intensity I in Figs. 1 (right part), 2, and 3. For a Raman shift ¯ h ω R , at light excitation energy ¯ h ω L , I is asfollows : I αβ ( ω R , ω L ) = ¯ h ( ω L − ω R ) ( π ) c N q ∑ m q " ∂ ε αβ ( ω L ) ∂ U m q ∂ U m ( − q ) n ( ω m q ) + ω m q δ ( ω R − ω m q ) ) . (1)The tensorial nature of I ( α and β referring to the polarization of incoming and outgoing light) will not play an importantrole in our discussion. Apart prefactors, the Raman intensity I is essentially obtained by summing over contributions fromall phonons, characterized by their momentum q and branch index m . Such contributions are products of the change of thedielectric response ε due to each phonon and its time-reversed homolog, by a function that depends only on the phononfrequency ω m q and a Dirac-delta, that expresses the conservation of energy. This formula is explained in more detail in theMethods section, which provides its derivation, including the “overtone” approximation, as well as the description of dielectricresponse calculations, with or without excitonic effects.In Eq. (1), the Dirac-delta function, integrand of the two-phonon density of states g DOS ( ω R ) = N q ∑ m q δ ( ω R − ω m q ) (2)is actually weighted by a function of the phonon mode and excitation energy (the term in square brackets in Eq. (1)). Onecan expect such weighting factor to be a reasonably smooth function of the wavevector q . By contrast, the Dirac-delta functionof Eq. 2 needs to be properly integrated throughout the Brillouin zone. It induces van Hove singularities, that dominate theoverall shape of I .The second-order Raman scattering intensity can thus be decomposed into contributions coming from different points ofthe Brillouin zone for phonons. The assignment of peaks or characteristic features to selected phonons for Si had been estab-lished already in the seventies . Our methodology now allows one to examine accurately the behaviour of the intensity ofeach feature as a function of the light excitation energy ¯ h ω L , especially close to the resonance, as first attempted in the ninetiesby Grein and Cardona using the theoretical and computational tools available at that time. This is shown in Fig. 5, wherethe evolution of the diagonal component of the IPA Raman spectrum with respect to the excitation energy is depicted for eightcharacteristic wavevectors in the Brillouin zone. In this framework, ‘optical’ transitions between different electronic wavevec-tors, i.e., between k and k + q with non-zero q are forbidden in the equilibrium structure. The phonon vibrations, however,may change some of these transitions from forbidden to allowed, giving rise to non-zero momentum matrix elements, thusleading to finite second-order derivatives of the dielectric function with respect to the phonon displacement. As a reference,the electronic and phononic band structures are shown in Fig. 6, along with the one-phonon density of states (related to thetwo-phonon density of state given in Eq. (2), by a simple scaling of the phonon frequency axis). As shown in Fig. 6 (b), at the X point, the highest phonon frequency is about 470 cm − , twice this value being 940 cm − ; the electronic resonance energyis 1.24 eV. At the L point, twice the highest phonon frequency is slightly larger around 1000 cm − , with a resonance at 2.14eV. Likewise, for the Γ point, the corresponding values are 1040 cm − and 3.2 eV. From a careful inspection of Fig. 5, wecan conclude that the peak structures in Fig. 1 are related to these resonances. To highlight their behavior as a function of theexcitation energy, we plot in Fig. 7 the second derivatives of the dielectric function at various points of the BZ. At the X point,it shows a clear maximum at the energy of the indirect band gap. The four curves related to the L point do not peak at theirresonance frequency, but increase with excitation energy, with a reduced rate above 2 eV. The curve at Γ grows unabated evenbeyond 3.0 eV. From such different characteristics, one can also conclude that the two-phonon density of states alone is notsufficient to describe second-order Raman spectra in general, and even more so in resonance condition.To assess the impact of excitonic effects, we have performed the corresponding calculations with the derivatives of thedielectric function for phonon wavevectors Γ , X , and L obtained by solving the BSE. A comparison with the IPA results isgiven in Fig. 8, where we emphasize the logarithmic scale of the vertical axis. The corresponding Raman spectra obtainedfor several excitation energies are presented in Fig. 3. As a general trend, the ratio between the BSE and IPA intensities(highlighted in the inset) increases with the excitation energy, exhibiting a maximum at the direct band gap. This effect is .00.10.20.30.40.50.60.70.8 × × × × − ) Raman shift (cm − ) I n t e n s i t y ( a r b . un i t s ) I n t e n s i t y ( a r b . un i t s ) I n t e n s i t y ( a r b . un i t s ) I n t e n s i t y ( a r b . un i t s ) Figure 5. (Color online). Diagonal component of the second-order Raman spectrum obtained in the IPA (blue lines), and itswavevector decomposition for four different excitation energies ¯ h ω L : 1.58 eV, 1.96 eV, 2.33 eV and 2.80 eV. Eight q wavevectors are considered, in reduced coordinates : Γ = ( , , ) , L = ( / , / , / ) , X = ( / , / , ) , W = ( / , / , / ) , (1) = ( / , / , ) , (2) = ( / , / , ) , (3) = ( − / , / , ) , (4) = ( / , / , / ) . The two-phonondensity of states is also represented for sake of comparison (red lines). E n e r g y ( e V ) ω ( c m − ) Figure 6. (Color online). (a) Electronic Kohn-Sham band structure of silicon, with a scissor shift of 0.65 eV for theconduction bands. Possible direct transitions at Γ and indirect transitions from Γ to X and from Γ to L are indicated by thegreen, red, and blue arrows, respectively. (b) Phonon band structure and (c) phonon density of states.similar to the first-order enhancement effect reported in Ref.34. A very strong excitonic effect can be observed between 2 and3 eV for the third and fourth mode of the L point, as depicted in Fig. 8 for the latter. It can be explained by the presence ofbound exciton states (inside the gap) with finite momentum that become visible at second order due to the presence of finite-momentum phonons. The ratio between BSE and IPA results dramatically increases in this energy range due to the changedonset of the second-order response. The other L -point modes, around 100 and 500 cm − , do not exhibit strong coupling withexcitons. Overall, exciton-phonon coupling results in a pronounced intensity increase in the part of the Raman spectrum whichis dominated by L modes, for excitation energies between 2 and 3 eV and Raman shifts between 600 and 900 cm − (Fig. 3).This result is a central outcome of the present work.For the other modes, indirect transitions appear less important for the resonance behavior than direct transitions. Forexample, in the regions around 1000 and 1200 cm − , all the intensities increase with energy with a similar rate. Since therelative intensities are well reproduced within the IPA, this explains why the IPA spectra in this range are in already quite goodagreement with experiment (see Fig. 1 for a comparison between IPA and BSE in this range). On the other hand, an accuratedescription of absolute intensities requires the inclusion of excitonic effects as they lead to changes of one order of magnitudefor excitation energies of around 3 eV, as already shown in Fig. 2.Silicon, a paradigmatic material in solid-state physics, might however not be the material in which exciton-phonon cou-pling has the biggest impact on the second-order Raman spectrum. Indeed, for materials that possess infra-red active phononmodes (unlike silicon), the Fr¨ohlich coupling with electrons diverges for small wavevectors, and dominates the exciton-phononinteraction. This has been explored for the second-order Raman case on the basis of model Hamiltonians in the eighties andnineties, see e.g. Ref. 24 and references therein. The present approach can be generalized to the treatment of such Fr¨ohlichcoupling. Conclusions
Thanks to a new method that combines frozen-phonon supercell calculations of dielectric functions and their derivatives withphonon spectra from density-functional perturbation theory, we have been able to compute frequency-dependent second-orderRaman spectra. We have applied our approach to the computation of the second-order Raman intensities of silicon for aseries of excitation energies. We not only reach excellent agreement with measured spectra but provide insight into theresonance behavior through analysis in terms of electronic structure and contributions from phonons at different points of theBrillouin zone. Excitonic effects turn out to be crucial to determine absolute Raman intensities as already found for the first-order scattering . Moreover, we have shown that strong excitonic effects may also appear selectively for some vibrationalmodes of different phonon wavevectors, modifying the relative intensities of different peaks. Our approach leads to a deeperunderstanding of the physics involved in second-order Raman processes and facilitates the interpretation of experimental .0 0.5 1.0 1.5 2.0 2.5 3.010 − − − − − − ΓL 1+2L 3L 4L 5+6X 1+2X 3+4X 5+6 1.24 eV 2.14 eVExcitation energy ~ ω L (eV) | ∂ ε/dU | Figure 7. (Color online). Second derivatives of the IPA dielectric function with respect to the excitation energy for differentmodes at the Γ , X , and L points. The notation “ q m” represents the m th branch of a q-point ( Γ , X , or L ), while “q m + n ”represents the sum of the corresponding derivatives of the degenerate branches m and n at this point. The electronic transitionthresholds for X (1.24 eV) and L (2.14 eV) are also indicated.results. While in this article we have focused on silicon—a simple but prototypical material—our methodology is generallyapplicable to a broad range of materials. Methods
We focus first on the theoretical methods, and then we give the computational details.In the harmonic approximation for the vibrational properties of periodic solids, the following expression can be used tocompute the Raman spectrum for two Stokes processes: I αβ ( ω R , ω L ) = ( ω L − ω R ) ¯ h ( π ) c N q ∑ q ∑ m , m " ∂ ε αβ ( ω L ) ∂ U m q ∂ U m ( − q ) ω m q ω m ( − q ) [ n m q + ][ n m ( − q ) + ] δ γ ( ω R − ( ω m q + ω m ( − q ) ) ) ) , (3)where α and β are Cartesian indices related to the polarization of the incoming and outgoing light, q are phonon wavevectorssampling the full BZ homogeneously (with a total of N q points), m and m are phonon branch indices, ω R is the spectralfrequency (Raman shift), ω L is the excitation frequency, c the speed of light, ω m q is the phonon frequency of the m th phononbranch at q , n m q is the temperature-dependent phonon occupation factor n m q = ( e ¯ h ω m q kT − ) − and δ γ is a Dirac-delta functionreplaced by a Lorentzian function of width γ for numerical purposes. In Eq. (3) , the derivative of the dielectric function, ε , is obtained through an expansion of the dielectric function with respect to the phonon eigendisplacements. Extracting thesecond-order process from this expansion and using the adiabatic approximation gives, for the electronic part, the followingsecond-order derivative: ∂ ε αβ ∂ U m q ∂ U m ( − q ) = ∑ κα l , κ ′ β l ′ ∂ ε αβ ∂ R l κα ∂ R l ′ κ ′ β e i q ( R l − R l ′ ) U m q κα U m ( − q ) κ ′ β , (4)where R l is the coordinate of the l th periodic cell, R l κα the coordinate of the κ th atom along α axis in the l th periodic cell. U m q κα are the atomic eigendisplacements of the dynamical problem, following the conventions of Ref. 43. -BSEΓ-IPA L4-BSEL4-IPA L5+6-BSEL5+6-IPA X5+6-BSEX5+6-IPAL4 X5+6L5+6Γ010 3.00.0 1.0 2.0 − − − − − − − − I BSE /I IPA | ∂ ε/dU | Excitation energy ~ ω L (eV) Figure 8. (Color online). Second derivatives of the dielectric function with respect to the excitation energy for thehighest-energy modes at Γ , X and L as well as the fourth mode at the L point. The inset shows the ratio between the BSE andthe IPA results.The numerical evaluation of Eq. 3 is highly challenging since it requires careful integration over the Brillouin zone andsummation over two mode indices. This contrasts with the summation over only one mode index, at the zone center, requiredfor first-order Raman intensities.We now describe how to numerically evaluate from first-principles the above-mentioned equations. A full ab-initio phononband structure can be obtained from DFPT . As a first step, the phonon frequencies ω m q and the phonon eigendisplace-ments U m q κα are computed for a set of homogeneous q -points, and interpolated on a denser mesh of wavevectors thanks tothe knowledge of interatomic force constants, as described in Ref. 43. Unfortunately, existing DFPT implementations do notprovide second derivatives of the frequency-dependent dielectric function with respect to atomic positions especially whenexcitonic effects are to be included. We thus numerically differentiate the dielectric function with respect to the eigendis-placements, which requires supercells of different size, depending on the phonon wavelength. In the frozen-phonon approach,a supercell, compatible with the selected wavevector q is generated: the axes of this supercell (cid:8) R ′ j (cid:9) fulfill the constraint e i qR ′ j =
1, which means that q is a reciprocal vector of the supercell. The atoms of the supercell are then displaced with realdisplacements defined from the phonon eigendisplacements and phase factors u κα l , (cid:2) U κα ( q ) e i qR l + U κα ( − q ) e − i qR l (cid:3) , v κα l , − i (cid:2) U κα ( q ) e i qR l − U κα ( − q ) e − i qR l (cid:3) . (5)In the supercell, k + q points of the Brillouin zone are folded back to k points of the supercell. This allows for probing indirecttransitions.The calculations of the optical spectrum can be performed with different levels of accuracy, ranging from the indepen-dent particle approximation to time-dependent DFT (TDDFT) and, finally, many-body perturbation theory (MBPT) with theinclusion of excitonic effects by means of the Bethe-Salpeter equation .We now detail why even for a small system like silicon such calculations are rather demanding (extremely demanding inthe BSE case). First, for any wavevector, the summations over the N modes phonon modes in Eq. 3 require N modes × N modes evaluations of the derivatives. If the dielectric tensor is evaluated on a grid of 5 × × = × =
30 evaluations of the dielectric function to reachthe final Raman spectrum of silicon. This number can be even reduced to only 5 by using symmetry arguments: acousticmodes do not contribute and the three other modes are degenerate.It should also be noted that, unlike first-order spectra, second-order resonant Raman scattering requires the use of super-cells commensurate with the wavevectors, and the computational cost increases quickly with the size of the supercell. Indeed,the number of bands included in the calculations scales linearly with the number of atoms. In the case of IPA, the numberof electronic transitions scales linearly with the product of the number of valence and conduction bands and linearly with thenumber of k-points. Taking into account the folding of bands in the case of supercell, and the computational cost required for he dipole matrix elements, the overall scaling is roughly cubic with the size of the supercell. In BSE, however, the kernelintroduces coupling between transitions, and the number of matrix elements in the transition basis set scales with the squareof the product of the number of valence and conduction bands, and with the square of the number of k-points. Because ofthe reduction of the number of k-points with larger supercells, the number of matrix elements therefore increases with thesquare of the size of the supercell. Since the evaluation of the screened Coulomb matrix elements requires a double sum overplanewaves, we conclude that the ratio between a BSE run for a supercell and the one for a unit cell increases as the fourthpower of the number of replicas in the supercell. Finally, as already observed in Ref. 34, the derivatives of the dielectricfunction are highly sensitive to the sampling of the BZ, leading to very demanding calculations in the case of BSE. Even ifthe frozen-phonon methodology requires a large amount of evaluations of dielectric functions, we have succeeded to applythe technique to three-dimensional materials by using approximations and numerical techniques as follows.As mentioned in the Results section, the so-called “overtone” contributions, i.e., considering only derivatives coming fromtwice the same phonon mode, except when degenerate, should largely dominate other second-order contributions. In Eq. (3)this corresponds to neglecting all terms where ω m = ω m and yields Eq. (1).The overtone approximation allowed us to reduce significantly the number of required computations. Eventually, only one-dimensional derivatives are needed for non-degenerate modes, reducing the number of evaluations required for one wavevectorfrom 900 for the full expression to 5 × =
30. For high-symmetry points of the BZ, few two-dimensional derivatives are stillrequired between two different modes belonging to the same degenerate subspace. As shown in Fig. 1(b), the frequency-dependent spectrum is well reproduced for silicon with this approximation. This figure corresponds to the Z ( XX ) Z scatteringconfiguration.In order to further reduce the computational load, the second-order derivatives of the dielectric function have been eval-uated only on a coarse mesh of phonons wavevectors, while the phonon frequencies, the occupation numbers and Diracfunctions have been calculated via DFPT techniques and summed on a dense mesh. This is equivalent to taking an averageddensity of states around the coarse points as the density function of Ref. 27.The corresponding grid parameters are now described, as well as other computational details. We use the ABINIT pack-age , with Troullier-Martin pseudopotentials (available from the ABINIT package ). Relaxed cell parameters of 10.20Bohr are used with cut-off values for the wavefunctions of 8 Hartree. The DFPT phonons displacements and frequenciesare interpolated from a coarse Γ -centered phonon grid of 8 × × × ×
70 mesh, to obtain the phonondensity of states, following Ref. 43, while a 4-times shifted 8 × × , the dielectric function is evaluated for different supercellscorresponding to a non-shifted 4 × × × ×
24 grid for the primitive cell is used. In the supercell case, the Brillouin zone is folded,and the integration grid adapted accordingly. The convergence has been checked on the final Raman spectrum with respect toa 4-times shifted 2 × × Γ , X and L wavevectors (corresponding to a 2 × × Γ -centered grid). Dielectricfunctions obtained on 12 × ×
12 wavevectors in the BZ have been averaged following the technique described in Ref.34 toreach samplings equivalent to 24 × ×
24. 3 valence bands and 6 conduction bands were used for the primitive cells, theirnumbers being increased proportionally with the number of atoms in the supercells. The model dielectric function describedin Ref.57 with ε ∞ =
12 has been used to calculate the screening matrix elements. A cut-off energy of 4 Ha is used to describethe screening matrix in reciprocal space.In order to present a Raman spectrum obtained with a sufficient number of phonon frequencies in the BSE framework,we have interpolated the ratio between BSE and IPA intensities towards a non-shifted 4 × × ( ω L − ω R ) factor is discarded. This factor is present also in the Ramanspectrum of reference wide-band gap materials, like CaF , to which the comparison is made in order to establish the absoluteintensity. Hence, it can be discarded. Other factors are expected to be constant in wide band-gap materials. References Yu, P. Y. & Cardona, M.
Fundamentals of Semiconductors: Physics and Material Properties (Springer, 2010). Weber, W. & Merlin, R.
Raman Scattering in Materials Science (Springer, 2000). Cardona, M. & Guntherodt, G.
Light scattering in solids / edited by M. Cardona ; with contributions by M.H. Brodsky ...[et al.] (Springer-Verlag Berlin ; New York, 1975). . Cardona, M. & Guntherodt, G.
Light scattering in solids II / edited by M. Cardona ; with contributions by M. Cardona ...[et al.] (Springer-Verlag Berlin ; New York, 1982). Lanzani, G.
The Photophysics behind Photovoltaics and Photonics (Wiley, 2012). Pan, Z. et al.
Polarization-resolved spectroscopy imaging of grain boundaries and optical excitations in crystalline organicthin films.
Nat. Commun. , 10.1038/ncomms9201 (2015). Akselrod, G. M. et al.
Visualization of exciton transport in ordered and disordered molecular solids.
Nat. Commun. ,10.1038/ncomms4646 (2014). Nozik, A. J. Photovoltaics: Separating multiple excitons.
Nat. Photon. , 272–273 (2012). Ferrari, A. C. & Basko, D. M. Raman spectroscopy as a versatile tool for studying the properties of graphene.
Nat. Nano. , 235–246 (2013). Thomsen, C. & Reich, S. Double resonant Raman scattering in graphite.
Phys. Rev. Lett. , 5214–5217 (2000). Berkdemir, A. et al.
Identification of individual and few layers of WS using Raman spectroscopy. Sci. Rep. , 1755(2013). del Corro, E. et al. Excited excitonic states in 1l, 2l, 3l, and bulk WSe observed by resonant Raman spectroscopy. ACSNano , 9629–9635 (2014). Gaur, A., Sahoo, S., Scott, J. & Katiyar, R. Electron–phonon interaction and double-resonance Raman studies in mono-layer WS . J. Phys. Chem. C , 5146–5151 (2015).
Wang, G. et al.
Double resonant Raman scattering and valley coherence generation in monolayer WSe . Phys. Rev. Lett. , 117401 (2015). del Corro, E. et al.
Atypical exciton-phonon interactions in WS and WSe monolayers revealed by resonance Ramanspectroscopy. Nano Lett. , 2363–2368 (2016). Venezuela, P., Lazzeri, M. & Mauri, F. Theory of double-resonant Raman spectra in graphene: Intensity and line shapeof defect-induced and two-phonon bands.
Phys. Rev. B , 035433 (2011). Herziger, F. et al.
Two-dimensional analysis of the double-resonant 2d Raman mode in bilayer graphene.
Phys. Rev. Lett. , 187401 (2014).
Ferrari, A. C. & Robertson, J. Raman spectroscopy in carbons: from nanotubes to diamond.
Phil. Trans. R. Soc. A ,2267–2565 (2004).
Weber, M. C. et al.
Temperature evolution of the band gap in BiFeO traced by resonant Raman scattering. Phys. Rev. B , 125204 (2016). Go, S., Bilz, H. & Cardona, M. Bond charge, bond polarizability, and phonon spectra in semiconductors.
Phys. Rev. Lett. , 580 (1975). Cardona, M. & Allen, P.B. Theory of the resonant Raman scattering by two phonons in germanium and silicon.
HelveticaPhysica Acta , 307 (1985). Trallero-Giner C., Cantarero A. & Cardona, M. One-phonon resonant Raman scattering: Fr¨olich exciton-phonon interac-tion.
Phys. Rev. B , 4030 (1989). Grein C.H., Zollner S. & Cardona, M. Microscopic theory of second-order Raman scattering in silicon under uniaxialstress..
Phys. Rev. B , 6633 (1991). Garcia-Cristobal, A., Cantarero, A., Trallero-Giner C. & Cardona, M. Excitonic model for second-order resonant Ramanscattering.
Phys. Rev. B , 13430 (1994). Strauch, D. et al.
Full ab initio calculation of second-order infrared and Raman spectra of elemental semiconductors.
Physica B , 442 (1996).
Renucci, J. B., Tyte, R. N. & Cardona, M. Resonant Raman scattering in silicon.
Phys. Rev. B , 3885–3895 (1975). Knoll, P. & Ambrosch-Draxl, C. Raman scattering of atomic vibrations in anharmonic potentials. In
Proceedings of theInternational Workshop on Anharmonic Properties of High - Tc Cuprates , 220 (World Scientific Publishing Co, Singapore,1995). Eds. Eds. D. Mihailovic, G. Ruani, E. Kalids, and K. A. M¨uller.
Ambrosch-Draxl, C. et al.
Raman scattering in YBa Cu O : a comprehensive theoretical study in comparison withexperiments. Phys. Rev. B , 064501 (2002). Puschnig, P., Ambrosch-Draxl, C., Henn, R. W. & Simon, A. Electronic properties and Raman spectra of rare-earthcarbide halides investigated from first principles.
Phys. Rev. B , 024519 (2001). Ravindran, P. et al.
Raman- and infrared-active phonons in superconducting and non-superconducting rare-earthtransition-metal borocarbides from full-potential calculations.
Phys. Rev. B , 104507 (2003). Spitaler, J., Sherman, E. Y. & Ambrosch-Draxl, C. Zone-center phonons in NaV O : A comprehensive ab-initio studyincluding Raman spectra and electron-phonon interaction. Phys. Rev. B , 014302 (2007). Spitaler, J., Sherman, E. Y. & Ambrosch-Draxl, C. First-principles study of phonons, optical properties, and Ramanspectra in MgV O . Phys. Rev. B , 064304 (2008). Spitaler, J., Sherman, E. Y., Ambrosch-Draxl, C. & Evertz, H. G. Lattice vibrations in CaV O and their manifestations:A theoretical study based on density functional theory. New J. Phys. 11 , 113009 (2009). Gillet, Y., Giantomassi, M. & Gonze, X. First-principles study of excitonic effects in raman intensities.
Phys. Rev. B ,094305 (2013). Baroni, S. & Resta, R.
Ab initio calculation of the low-frequency raman cross section in silicon.
Phys. Rev. B ,5969–5971 (1986). Baroni, S., Giannozzi, P. & Testa, A. Green’s-function approach to linear response in solids.
Phys. Rev. Lett. , 1861–1864 (1987). Gonze, X. Adiabatic density-functional perturbation theory.
Phys. Rev. A , 1096–1114 (1995). Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P. Phonons and related crystal properties from density-functionalperturbation theory.
Rev. Mod. Phys. , 515–562 (2001). Veithen, M., Gonze, X. & Ghosez, P. First-principles study of the electro-optic effect in ferroelectric oxides.
Phys. Rev.Lett. , 187401 (2004). Veithen, M., Gonze, X. & Ghosez, P. Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors fromfirst-principles density functional perturbation theory.
Phys. Rev. B , 125107 (2005). Gonze, X., Rignanese, G.-M. & Caracas, R. First-principle studies of the lattice dynamics of crystals, and related proper-ties.
Z. Kristallogr. , 458 (2005).
Caracas, R. & Bobocioiu, E. The wurm project—a freely available web-based repository of computed physical data forminerals.
Am. Mineral. , 437 (2011). Gonze, X. & Lee, C. Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic forceconstants from density-functional perturbation theory.
Phys. Rev. B , 10355–10368 (1997). Temple, P.A. & Hathaway, C.E. Multiphonon Raman spectrum of silicon.
Phys. Rev. B , 3685 (1973). Wang, C.S., Chen, J.M., Becker, R. & Zdetsis A. Second order Raman spectrum and phonon density of states of silicon.
Phys. Lett. , 517 (1973).
Trommer, R. & Cardona, M. Resonant Raman scattering by 2TO phonons and the ordering of conduction band minimain GaAs.
Solid State Comm. , 153 (1977). Pasternak, A., Cohen, E., & Gilat, G. Calculation of second-order Raman scattering for KBr, NaCl, and MgO crystals.
Phys. Rev. B , 551 (1974). Reich, S. et al.
Resonant Raman scattering in cubic and hexagonal boron nitride.
Phys. Rev. B , 205201 (1995). Gatti, M. & Sottile, F. Exciton dispersion from first principles.
Phys. Rev. B , 155113 (2013). Verdi, C. & Giustino, F. Fr¨ohlich Electron-Phonon Vertex from First Principles.
Phys. Rev. Lett. , 176401 (2015).
Hanke, W. & Sham, L. J. Many-particle effects in the optical excitations of a semiconductor.
Phys. Rev. Lett. , 387–390(1979). Onida, G., Reining, L. & Rubio, A. Electronic excitations: density-functional versus many-body green’s-function ap-proaches.
Rev. Mod. Phys. , 601–659 (2002). Gonze, X. et al.
Abinit : First-principles approach of materials and nanosystem properties.
Comput. Phys. Comm. ,2582 (2009).
Gonze, X. et al.
Recent developments in the ABINIT software package.
Comput. Phys. Comm. , 106 (2016).
Gonze, X. & Beuken, J.-M. ABINIT – abinit. URL . (2016) (Date of access: 02/10/2016). Ambrosch-Draxl, C. and Sofo, J. 0. Linear optical properties of solids within the full-potential linearized augmentedplanewave method.
Comp. Phys. Commun. , 1 (2006).
Cappellini, G., Del Sole, R., Reining, L. & Bechstedt, F. Model dielectric function for semiconductors.
Phys. Rev. B ,9892–9895 (1993). Acknowledgments
Y.G. and M.G. acknowledge financial support by the Fonds National de la Recherche Scientifique (FNRS, Belgium) and theESF through the ESF/Psi-k2 programme. C.D. and S.K. appreciate support from the German Science Foundation (DFG)through CRC 658. We thank Yann Pouillon and Jean-Michel Beuken for their valuable technical support and help withthe test and the build system of ABINIT. Computational resources have been provided by the supercomputing facilities ofthe Universit´e catholique de Louvain (CISM/UCL) and the Consortium des Equipements de Calcul Intensif en F´ed´erationWallonie Bruxelles (CECI) funded by the Fonds de la Recherche Scientifique de Belgique (FRS-FNRS). The present researchbenefited from computational resources made available on the Tier-1 supercomputer of the F´ed´eration Wallonie-Bruxelles,infrastructure funded by the Walloon Region under the grant agreement n o Author contributions
Y.G. and S.K. developed the methodology. Y.G. performed the calculations. Y.G., C. D., X. G. analyzed the data. S.K., M.G.,C.D., X.G. supervized the work. Y.G., M.G., C. D., X. G. wrote the paper.
Competing financial interests