Excitons on a microscopic level: The mixed dynamic structure factor

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Excitons on a microscopic level: The mixed dynamic structure factor Igor Reshetnyak, Matteo Gatti, Francesco Sottile, Lucia Reining

The dynamic structure factor of materials is proportional to their linear electronic response and it displays their excitation spectra. Usually the response is measured on the same length scale as the perturbation. Here we illustrate that much can be gained by studying also the mixed dynamic structure factor, which connects different spatial components of perturbation and response. We extend state-of-the-art ab initio calculations to access the mixed dynamic structure factor, including excitonic effects. Using bulk silicon and lithium fluoride as prototype examples, we show that these effects play a crucial role above and below the quasi-particle gap, and are needed in order to explain coherent inelastic x-ray scattering experiments. Our approach also yields important information concerning the microscopic structure of time-dependent density functional theory.
One of the key concepts in condensed matter theory is screening, the modification of a potential felt by a charge due to the rearrangement of other charges [1]. In many cases the dominant contribution to this phenomenon can be described in linear response [2]. This means that the knowledge of the response of a system to an external perturbation is entirely contained in the density-density response function χ,which is determined only by the material.The poles of χ in frequency space are the excitation energies of the system. The densitydensity response function is measured directly or indirectly in many spectroscopy experiments [3], such as electron energy loss spectroscopy (EELS) [4], optical absorption [5,6], or inelastic x-ray scattering (IXS) [7], which yields the dynamic structure factor (DSF) that is proportional to the imaginary part of χ. Knowing χ can also help to understand or predict effects such as a significant local rearrangement of charges as response of the system to a perturbation, which can have dramatic effects on the structure, for example self-trapped excitons [8]. Finally, screening is also one of the fundamental processes that govern the behavior of all many-body systems. Therefore it appears naturally as a building block in the formulation of many-body perturbation theory [9], via the screened Coulomb interaction W = v c + v c χv c , where v c is the bare Coulomb interaction. The widely-used GW approximation (GWA) [10], for example, uses W as effective interaction.
Optics, EELS and IXS probe the response of the system on the same length scale as the perturbation. Besides the spectroscopic information, this allows one also to reconstruct charge excitations as an averaged function of space and time from experimental IXS spectra [51][52][53][54][55]. However, in an inhomogeneous material even a spatially monochromatic perturbation creates a response on different length scales [56,57], depending on the local structure of the material. Only when all these components of the response are known can one describe induced charges with spatial resolution, and can one determine important many-body effects that depend on all the components of W . For these reasons, it is highly desirable to extend state-of-the-art advanced theoretical and numerical approaches to the description of the full inhomogeneous response of materials including excitonic effects. This gives access to the mixed DSF (MDSF), which is a matrix in reciprocal space whose diagonal is the ordinary DSF [7]. The MDSF can be measured by coherent IXS (CIXS) [7]. On the theory side, some off-diagonal elements of the MDSF in silicon were calculated in the adiabatic local density approximation (ALDA) of time-dependent den-sity functional theory (TDDFT) [58,59], which is similar to the RPA, whereas to the best of our knowledge excitonic effects have never been accessed.
The aim of the present work is to extend ab initio BSE calculations to access the MDSF including excitonic effects, to benchmark against experiment, analyze and make predictions. This also leads to new guidelines for modeling approximations to TDDFT.
In periodic crystals the density-density response function can be written in reciprocal space as a function and matrix χ GG (q, ω), where q belongs to the first Brillouin zone and G, G are reciprocal-lattice vectors. It yields the induced charge as response to an external potential In a homogeneous material χ is diagonal. In an inhomogeneous material the diagonal of χ yields only a spatially averaged response, and not its local values.
To access excitonic effects, we have to extend the BSE approach to the calculation of all elements χ GG (q, ω) of the response matrix. In its state-of-the-art ab initio version this approach uses a quasiparticle (QP) approximation, and replaces the electron-hole interaction by a statically screened Coulomb interaction W (ω = 0). This allows one to express χ in terms of the eigenvectors A λ and eigenvalues E λ of the two-particle hamiltonian H exc =Ĥ el +Ĥ h +v −Ŵ , whereĤ el andĤ h contain the occupied and empty QP bands calculated in the GWA,v is the electron-hole exchange interaction due to the bare v c , and −Ŵ is the screened direct electron-hole attraction that is responsible for excitonic effects [26]. Settinĝ W = 0 corresponds to performing an RPA calculation that uses as input the GW bandstructure: we will refer to this as RPA+GW. Using instead the LDA bandstructure as input leads to the RPA+LDA approximation, which is often simply called RPA. For a system with a gap, one can express the full matrix χ in terms of the eigenvalues E λ and eigenvectors A λ ofĤ exc as: where O λλ is the overlap matrix of the coefficients A λ that mix the transitions t among QP states (n → n bands and k − q → k wavevectors), with oscillator strengthρ t (q + G) = φ nk−q | e −i(q+G)·r |φ n k (for more details see Ref. [45], where the diagonal G = G of Eq.(1) was calculated).
to determine the plasmon band gap opening at the Brillouin zone boundary due to the periodic lattice potential [60,[69][70][71][72]. This allows us to benchmark our theoretical development against experiment for silicon, before making a prediction for the case of a strongly bound exciton in LiF.
In both these centrosymmetric crystals (see [73] for the general formulation) the MDSF S GG (q, ω) is real valued and related to χ via: S GG (q, ω) = −[Imχ GG (q, ω) + Imχ G G (q, ω)]/(2π). Fig. 1(a) shows the calculated S GG (q, ω)of Si at q = (−1/2, −1/2, −1/2), G = (0, 0, 0) with either the diagonal element G = G or the offdiagonal one with G = (1, 1, 1),compared to the experimental results from [22,60,61] [74]. The diagonal element is characterised by a broad plasmon peak. As it is found in general [39], the RPA+GW calculation overestimates the plasmon energy, and the electron-hole attraction yields a shift to lower energies, close to 20 eV, in good agreement with experiments [22,61]. The off-diagonal elements of S GG (q, ω) look qualitatively different from the diagonal ones. In particular, their spectrum can even be negative [58,59]. Our results show that many-body effects act in a similar way as in the diagonal elements: compared to experiment, the RPA+GW gives a qualitatively correct description, but it overestimates the energy of all structures seen in experiment. The electron-hole attraction in the BSE shifts spectral weight to lower energies, so that there is a first peak at 13.8 eV and a dip at 19.6 eV. Moreover, the first peak is sharper in BSE than in RPA+GW. This brings both the spectral shape and the position of the peaks much closer to experiment [60] than the RPA+GW results. One may find it intuitive that the BSE acts in a similar way on diagonal and off-diagonal elements. However, this is far from trivial: excitonic effects in silicon are interference effects, based on a subtle interplay of intensities and phases, not simply a shift of energies [36,75]. Without knowing the result, one could not have predicted the picture of Fig. 1(a).
The very good agreement between theory and experiment for silicon, and the significant excitonic effects on off-diagonal elements that emerge from the results, give strong motivation to explore and predict spectra for a material with much stronger excitonic effects. Fig. 1(b) shows the calculated MDSF of LiF for the same choice of q, G, G as above for silicon. Our calculated diagonal element agrees with literature results [42][43][44][45]; in particular, a strongly bound exciton at 14.6 eV, well within the QP gap given by the onset of RPA+GW, is observed in the measured spectra and correctly described by the BSE. The question is how the exciton will affect the offdiagonal elements. Our calculation shows that, as in the case of the diagonal, RPA+GW yields a spectrum that is quite featureless. It is negative and of significant intensity, but more than a factor of two weaker than the diagonal one. Including now the electron-hole interaction, the structures in the off-diagonal element above 15 eV reflect structures in the diagonal element with positive or negative sign and, as in the case of RPA+GW, clearly reduced intensity. Below the GW gap also in the offdiagonal element a strongly bound exciton peak appears. Strikingly, it is approximately the specular negative of its diagonal counterpart at 14.6 eV, with a very similar intensity. This shows that it could be detrimental to neglect off-diagonal elements, for example to get insight into bound excitons in real space [53].
The CIXS experiment is difficult, and published results are today still limited to the first pioneering works on silicon [60]. Our results, while showing the importance of an ab initio theoretical approach with predictive power, strongly motivate a renewed experimental effort. At the same time, they can impact other fields of theory. Indeed, in principle the full χ including excitonic effects can also be calculated in TDDFT [76,77]. With respect to the two-particle BSE, TDDFT has the computational advantage of relying on the oneparticle Kohn-Sham equation. It is therefore a widely used approach all over physics and chemistry [78][79][80][81]. In TDDFT, χ can be obtained from the linear response ma-trix equation χ = χ RP A + χ RP A f xc χ, where χ RP A is the RPA response function based on Kohn-Sham ingredients, and the exchange-correlation kernel f xc is the functional derivative of the Kohn-Sham exchange-correlation potential with respect to the density [26]. Like χ GG (q, ω), in a crystal χ RP A and f xc are matrices in reciprocal space. Unfortunately, the widely used approximations for f xc , such as the ALDA, fail to describe excitonic effects in extended systems [82]. More advanced kernels have been derived from the BSE [83][84][85][86][87][88][89][90][91][92] or modeled using related knowledge [93][94][95][96][97][98][99], but no known expression is to date able to well describe spectra including bound and continuum excitons with a reasonable computational effort.
One obstacle for approximating f xc is the fact that little is known, in particular concerning off-diagonal elements. If χ was given, in principle f xc could be studied by inverting the linear response equation, f xcGG (q, ω) = ω). However, up to now offdiagonal elements of χ GG (q, ω) in the presence of excitonic effects were unknown, and as Fig. 1 shows, they cannot be neglected.
The present work overcomes this difficulty, since we now have the full matrix χ GG (q, ω) and can therefore invert the linear response equation. In line with literature, we will distinguish two definitions: the f xc from TDDFT as explained above, and f mb ω). The latter "manybody" kernel does not have the difficult task to open the quasi-particle gap, because χ RP A+GW is built with the GW band structure. It is therefore a common starting point for successful model kernels (all kernels in Refs. [83][84][85][86][87][88][89][90][93][94][95][96][97] use this approach). In particular, we can address the following questions: What is the global structure of these kernels? Are there significant differences that might make one or the other easier to model? Can one detect a special feature linked to the bound exciton? Fig. 2(a) shows for LiF the real part of the q → 0 head elements G=G =0 of f xc and f mb xc . Both kernels have a pronounced frequency dependence. For this reason, the real parts must be accompanied by a respective imaginary part, shown in Fig. 2(b). For both f xc and f mb xc the imaginary parts are zero at low frequency. The onset of intensity is above 8 eV for f xc , but only above 12 eV for f mb xc . The reason lies in the fact that Im f xc has to destroy intensity in the absorption spectrum between the onset of the Kohn-Sham spectrum and the bound exciton, whereas f mb xc only has to create the bound exciton in the QP gap that is already open in χ RP A+GW [100]. The two kernels f xc and f mb xc start with opposite sign, both in the real and in the imaginary part, reflecting the fact that their task is to transfer oscillator strength to higher or lower energy, respectively.
Our calculations allow us to estimate the importance of the matrix components of the kernels beyond the only head. To this aim, Fig. 3 shows absorption spectra eval- uated from Im (ω) = Im 1 + v c χ(q → 0, ω) −1 G=G =0 , where χ GG (q, ω) is calculated taking into account either the full matrix f xc or f mb xc (which both, by definition and also in practice, yield the reference BSE result and are therefore not shown), or only the respective head. The reference spectrum is the black line. Using only the head of the matrix f xc leads to a disaster, with a spectrum that is washed out by oscillations, with even negative spectral weight, and the complete absence of a bound exciton. Restricting the matrix to the only head element in f mb xc , instead, has a less dramatic effect: in particular, the spectral weight remains mainly positive, and the strongly bound exciton appears and is reasonably well positioned. This is a clear indication for the fact that it is much easier to model f mb xc than f xc . Still, the excitonic transfer of spectral weight to lower energies is severely overestimated, suppressing the spectrum at higher energies, similar to results of model scalar and static kernels [94,101] [102].
Here, instead, we explore a strongly bound exciton, with the example of LiF shown in Fig. 2. The effective kernels show similar structures as the head elements of their matrix version. However, the amplitude of both real and imaginary part [104] is smaller for the effective f ef f xc than for the head of f xc . This means that off-diagonal elements in the full response equation compensate to some extent. Altogether, our results suggest to direct efforts towards the modeling of mb or of effective scalar, rather than full, kernels for spectroscopy.
In conclusion, we have shown that the off-diagonal elements of the mixed dynamic structure factor carry important excitonic effects below and above the quasi-particle gap. This insight could be gained thanks to a generalization of the state-of-the-art ab initio solution of the Bethe-Salpeter equation. Our calculation of the continuum excitonic effects in silicon bring the spectra of off-diagonal elements in good agreement with experiments. Moreover, we predict that the strongly bound exciton in LiF that is seen in optics, EELS and IXS should also be observable in some off-diagonal elements of the MDSF. These offdiagonal elements have an intensity that is comparable to the intensity of diagonal elements to which experiments are usually limited, and cannot be neglected when one wants to access the full spatially resolved density response. This motivates a renewed experimental effort, especially concerning coherent IXS. Having access to the full response matrices also gives new theoretical impulses. In particular, it allowed us to derive information concern-ing various flavors of effective exchange-correlation kernels in the framework of time-dependent density functional theory, suggesting the ones most amenable to simple approximations. This research was supported by the European Research Council (ERC Grant Agreement n. 320971) and by a Marie Curie FP7 Integration Grant within the 7th European Union Framework Programme. Computational time was granted by GENCI (Project No. 544).