Ab Initio Approach to Second-order Resonant Raman Scattering Including Exciton-Phonon Interaction

Raman spectra obtained by the inelastic scattering of light by crystalline solids contain contributions from first-order vibrational processes (e.g. the emission or absorption of one phonon, a quantum of vibration) as well as higher-order processes with at least two phonons being involved. At second order, coupling with the entire phonon spectrum induces a response that may strongly 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 in resonance 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 in excellent agreement with experiment, to different vibrations and regions of the Brillouin zone. We show that exciton-phonon coupling, computed from first principles, indeed strongly affects the spectrum in resonance conditions. The ability to analyze second-order Raman spectra thus provides direct insight into this interaction.


Results
We first describe the normalized intensities for Raman shifts in the 900-1100 cm −1 range. Such results are presented in the right panel of Fig. 1. On both levels of theory, IPA and BSE, the evolution of the intensity of the Raman shift with increasing light frequency (or photon excitation energy) follows very closely the experimental data. While most of the spectral weight is localized around 950 cm −1 for the lowest excitation energy of 1.17 eV, the shoulder below 1000 cm −1 is getting more pronounced with increasing excitation energy and starts to dominate above 2.4 eV. A third peak around 1030 cm −1 that is hardly visible in the bottom curves, is particularly pronounced only for the highest excitation energy of 3.41 eV. In excellent agreement with experiment, its intensity is nearly the same as the peak at 950 cm −1 , and somewhat lower than the one at 980 cm −1 .
Second, we examine the integrated intensities. Figure 2 shows a comparison between the IPA and the BSE absolute intensity, divided by ω ω − ( ) L R 4 (see Methods section for more detail), integrated over the range between 900 and 1100 cm −1 , as a function of the excitation energy. Even if in this range of frequencies the relative intensities computed in the two formalisms are very similar as shown in Fig. 1, their absolute intensities differ by a factor 5 for excitation energies around 3 eV, as exhibited 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 absolute intensity, similar to the one shown in Fig. 2, Fig. 3 shows that the use of the BSE, that is, the inclusion of exciton-phonon coupling, results in a pronounced intensity increase in the part of the Raman spectrum between 600 and 900 cm −1 , especially for 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 is dominated by the emission of two phonons that are connected by time-reversal symmetry: same frequency, and time-reversed collective eigendisplacement and momentum. From a group-theoretical point of view, derivatives with respect to such phonon pairs will always be active and contribute to the fully symmetric component of the spectrum 3 . However, group-theoretical arguments do not insure the vanishing of other contributions for arbitrary phonon wavevectors. While this approximation is justified 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 derivatives of the dielectric constant obtained within Density-Functional Perturbation Theory (DFPT), where the electric field is used as a perturbation 43 . In this approach, the excitation energy vanishes, but unlike in the IPA, the local-field effects are taken into account, which makes its level of sophistication intermediate between IPA and BSE. Figure 4 shows the dominance of the derivatives by the overtone pair of modes (by around one order of magnitude with respect prefactor, see Eq. 1 has not been included, and the spectrum has been renormalized with respect to the value at ω = . 1 0 L  eV. Note the logarithmic scale of the intensity. The inset shows the ratio between integrated intensities obtained by the IPA and the BSE. to other contributions in the worst case), which fully justifies our approach. This 'overtone' approximation was a crucial step in making the BSE computation tractable, as discussed in the Methods section, and to be able to treat excitonic effects. The confirmation of the validity of the 'overtone' approximation for silicon, discussed in prior works, also for other simple semiconductors 4,20,21,[23][24][25][26][44][45][46] is thus another significant result of our work. However, there are several cases in which features can be assigned to combination modes, like in refs 16, 47 and 48. In such case, the identification of relevant combinations might be done in the static limit using DFPT, followed by the computation of the contributions from specific mode combinations in the relevant energy spectrum (as is 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 the underlying physics, especially the role of excitonic effects. For this purpose, we first present the expression that yields the theoretical Raman intensity I in Fig. 1 (right part), 2, and 3. For a Raman shift ω R , at light excitation energy ω L , I is as follows: The tensorial nature of I (α and β referring to the polarization of incoming and outgoing light) will not play an important role in our discussion. Apart prefactors, the Raman intensity I is essentially obtained by summing over contributions from all phonons, characterized by their momentum q and branch index m. Such contributions are products of the change of the dielectric response ε due to each phonon and its time-reversed homolog, by a function that depends only on the phonon frequency ω mq and a Dirac-delta, that expresses the conservation of energy. This formula is explained in more detail in the Methods section, which provides its derivation, including the "overtone" approximation, as well as the description of dielectric response calculations, with or without excitonic effects.
In Eq. (1), the Dirac-delta function, integrand of the two-phonon density of states is actually weighted by a function of the phonon mode and excitation energy (the term in square brackets in Eq. (1)) 2, 4 . One can expect such weighting factor to be a reasonably smooth function of the wavevector q. By contrast, the Dirac-delta function of Eq. 2 needs to be properly integrated throughout the Brillouin zone. It induces van Hove singularities, that dominate the overall shape of I.
The second-order Raman scattering intensity can thus be decomposed into contributions coming from different points of the Brillouin zone for phonons. The assignment of peaks or characteristic features to selected phonons for Si had been established already in the seventies 44,45 . Our methodology now allows one to examine accurately the behaviour of the intensity of each feature as a function of the light excitation energy ω L , especially close to the resonance, as first attempted in the nineties by Grein and Cardona 23 using the theoretical and computational tools available at that time. This is shown in Fig. 5, where the evolution of the diagonal component of the IPA Raman spectrum with respect to the excitation energy is depicted for eight characteristic wavevectors in the Brillouin zone. In this framework, 'optical' transitions between different electronic wavevectors, i.e., between k and k + q with non-zero q are forbidden in the equilibrium structure. The phonon vibrations, however, may ) obtained within density-functional perturbation theory. "Full" means that the emission of all possible phonon pairs are considered in the calculation, while for the "overtone" graph, only the emission of time-reversed phonon pairs are taken into account. change some of these transitions from forbidden to allowed, giving rise to non-zero momentum matrix elements, thus leading 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 the two-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 −1 , twice this value being 940 cm −1 ; the electronic resonance energy is 1.24 eV. At the L point, twice the highest phonon frequency is slightly larger around 1000 cm −1 , with a resonance at 2.14 eV. Likewise, for the Γ point, the corresponding values are 1040 cm −1 and 3.2 eV. From a careful inspection of Fig. 5, we can conclude that the peak structures in Fig. 1 are related to these resonances. To highlight their behavior as a function of the excitation 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 their resonance frequency, but increase with excitation energy, with a reduced rate above 2 eV. The curve at Γ grows unabated even beyond 3.0 eV. From such different characteristics, one can also conclude that the two-phonon density of states alone is not sufficient 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 the dielectric function for phonon wavevectors Γ, X, and L obtained by solving the BSE. A comparison with the IPA results is given in Fig. 8, where we emphasize the logarithmic scale of the vertical axis. The corresponding Raman spectra obtained for 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 similar to the first-order enhancement effect reported in ref. 34. A very strong excitonic effect can be observed between 2 and 3 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 of bound exciton states (inside the gap) with finite momentum 49 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 changed onset of the second-order response. The other L-point modes, around 100 and 500 cm −1 , do not exhibit strong coupling with excitons. Overall, exciton-phonon coupling results in a pronounced intensity increase in the part of the Raman spectrum which is dominated by L modes, for excitation energies between 2 and 3 eV and Raman shifts between 600 and 900 cm −1 (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. For example, in the regions around 1000 and 1200 cm −1 , all the intensities increase with energy with a similar rate. Since the relative intensities are well reproduced within the IPA, this explains why the IPA spectra in this range are in already quite good agreement with experiment (see Fig. 1 for a comparison between IPA and BSE in this range). On the other hand, an accurate description of absolute intensities requires the inclusion of excitonic effects as they lead to changes of one order of magnitude for 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 coupling has the biggest impact on the second-order Raman spectrum. Indeed, for materials that possess infra-red active phonon modes (unlike silicon), the Fröhlich coupling with electrons diverges for small wavevectors, and dominates the exciton-phonon interaction. This has been explored for the second-order Raman case on the basis of model Hamiltonians in the eighties and nineties, see e.g. ref. 24 and references therein. The present approach can be generalized to the treatment of such Fröhlich coupling 50 .

Conclusions
Thanks to a new method that combines frozen-phonon supercell calculations of dielectric functions and their derivatives with phonon spectra from density-functional perturbation theory, we have been able to compute frequency-dependent second-order Raman spectra. We have applied our approach to the computation of the  Second derivatives of the IPA dielectric function with respect to the excitation energy for different modes 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 transition thresholds for X (1.24 eV) and L (2.14 eV) are also indicated.
second-order Raman intensities of silicon for a series of excitation energies. We not only reach excellent agreement with measured spectra 26 but provide insight into the resonance behavior through analysis in terms of electronic structure and contributions from phonons at different points of the Brillouin zone. Excitonic effects turn out to be crucial to determine absolute Raman intensities as already found for the first-order scattering 34 . Moreover, we have shown that strong excitonic effects may also appear selectively for some vibrational modes of different phonon wavevectors, modifying the relative intensities of different peaks. Our approach leads to a deeper understanding of the physics involved in second-order Raman processes and facilitates the interpretation of experimental results. While in this article we have focused on silicon-a simple but prototypical material-our methodology is generally applicable 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 to compute the Raman spectrum 21,27 for two Stokes processes: where α and β are Cartesian indices related to the polarization of the incoming and outgoing light, q are phonon wavevectors sampling the full BZ homogeneously (with a total of N q points), m 1 and m 2 are phonon branch indices, ω R is the spectral frequency (Raman shift), ω L is the excitation frequency, c the speed of light, ω mq is the phonon frequency of the m th phonon branch at q, n mq is the temperature-dependent phonon occupation factor and δ γ is a Dirac-delta function replaced 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 the second-order process from this expansion and using the adiabatic approximation 27 gives, for the electronic part, the following second-order derivative: 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 mq are the atomic eigendisplacements of the dynamical problem, following the conventions of ref. 43. The numerical evaluation of Eq. 3 is highly challenging since it requires careful integration over the Brillouin zone and summation over two mode indices. This contrasts with the summation over only one mode index, at the zone center, required for first-order Raman intensities.
We now describe how to numerically evaluate from first-principles the above-mentioned equations. A full ab-initio phonon band structure can be obtained from DFPT 38,41 . As a first step, the phonon frequencies ω mq and the phonon eigendisplacements κα U mq are computed for a set of homogeneous q-points, and interpolated on a denser mesh of wavevectors thanks to the knowledge of interatomic force constants, as described in ref. 43. Unfortunately, existing DFPT implementations do not provide second derivatives of the frequency-dependent dielectric function with respect to atomic positions especially when excitonic effects are to be included. We thus numerically differentiate the dielectric function with respect to the eigendisplacements, 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 ′ R { } j fulfill the constraint ′ = e 1 iqR j , which means that q is a reciprocal vector of the supercell. The atoms of the supercell are then displaced with real displacements defined from the phonon eigendisplacements and phase factors In the supercell, k + q points of the Brillouin zone are folded back to k points of the supercell. This allows for probing indirect transitions.
The calculations of the optical spectrum can be performed with different levels of accuracy, ranging from the independent particle approximation to time-dependent DFT (TDDFT) and, finally, many-body perturbation theory (MBPT) with the inclusion of excitonic effects by means of the Bethe-Salpeter equation 51,52 .
We now detail why even for a small system like silicon such calculations are rather demanding (extremely demanding in the BSE case). First, for any wavevector, the summations over the N modes phonon modes in Eq. 3 require × N N modes m odes evaluations of the derivatives. If the dielectric tensor is evaluated on a grid of 5 × 5 pairs of displacements, in order to get the second-order derivatives, the number of calculations required for a single wavevector in the BZ becomes 25 × 6 2 = 900 evaluations of the dielectric tensor for silicon. In contrast, for first-order calculations, the number of computations for the zone-center wavevector amounts to 5 evaluations for each mode, thus 5 × 6 = 30 evaluations of the dielectric function to reach the final Raman spectrum of silicon. This number can be even reduced to only 5 by using symmetry arguments: acoustic modes 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 supercells 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 number of electronic transitions scales linearly with the product of the number of valence and conduction bands and linearly with the number of k-points. Taking into account the folding of bands in the case of supercell, and the computational cost required for the dipole matrix elements, the overall scaling is roughly cubic with the size of the supercell. In BSE, however, the kernel introduces coupling between transitions, and the number of matrix elements in the transition basis set scales with the square of the product of the number of valence and conduction bands, and with the square of the number of k-points. Because of the reduction of the number of k-points with larger supercells, the number of matrix elements therefore increases with the square of the size of the supercell. Since the evaluation of the screened Coulomb matrix elements requires a double sum over planewaves, we conclude that the ratio between a BSE run for a supercell and the one for a unit cell increases as the fourth power of the number of replicas in the supercell. Finally, as already observed in ref. 34, the derivatives of the dielectric function are highly sensitive to the sampling of the BZ, leading to very demanding calculations in the case of BSE. Even if the frozen-phonon methodology requires a large amount of evaluations of dielectric functions, we have succeeded to apply the 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 from twice 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 1 2 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 wavevector from 900 for the full expression to 5 × 6 = 30. For high-symmetry points of the BZ, few two-dimensional derivatives are still required 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 ( ) scattering configuration. In order to further reduce the computational load, the second-order derivatives of the dielectric function have been evaluated only on a coarse mesh of phonons wavevectors, while the phonon frequencies, the occupation numbers and Dirac functions have been calculated via DFPT techniques and summed on a dense mesh. This is equivalent to taking an averaged density 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 package 53,54 , with Troullier-Martin pseudopotentials (available from the ABINIT package 55 ). Relaxed cell parameters of 10.20 Bohr are used with cut-off values for the wavefunctions of 8 Hartree. The DFPT phonons displacements and frequencies are interpolated from a coarse Γ-centered phonon grid of 8 × 8 × 8 to a non-shifted 70 × 70 × 70 mesh, to obtain the phonon density of states, following ref. 43, while a 4-times shifted 8 × 8 × 8 sampling is used for the electrons.
When simulated in the IPA (without local-fields effects) 56 , the dielectric function is evaluated for different supercells corresponding to a non-shifted 4 × 4 × 4 phonon wavevector grid (containing 64 wavevectors), while for the electronic states, the equivalent of a centered 24 × 24 × 25 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 to a 4-times shifted 2 × 2 × 2 centered phonon wavevector grid (containing 32 wavevectors). A scissors operator of 0.65 eV is applied on top of the Kohn-Sham eigenvalues to mimic the opening of the gap by GW.
In order to assess the importance of excitonic effects, we have performed BSE calculations of the derivatives of the dielectric function for phonons with Γ, X and L wavevectors (corresponding to a 2 × 2 × 2 Γ-centered grid). Dielectric functions obtained on 12 × 12 × 12 wavevectors in the BZ have been averaged following the technique described in ref. 34 to reach samplings equivalent to 24 × 24 × 24. 3 valence bands and 6 conduction bands were used for the primitive cells, their numbers being increased proportionally with the number of atoms in the supercells. The model dielectric function described in ref. 57 with ε = ∞ 12 has been used to calculate the screening matrix elements. A cut-off energy of 4 Ha is used to describe the 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 × 4 × 4 k-point grid with a trilinear interpolation in reduced coordinates and used these ratios to correct the IPA intensities.
When integrated intensities are presented, the ω ω − ( ) L R 4 factor is discarded. This factor is present also in the Raman spectrum of reference wide-band gap materials, like CaF 2 , to which the comparison is made in order to establish the absolute intensity. Hence, it can be discarded. Other factors are expected to be constant in wide band-gap materials.