Correlation satellites in optical and loss spectra

Coupling of excitations leads to intriguing effects on the spectra of materials. We propose a cumulant formulation for neutral electronic excitations which opens the way to describe effects such as double plasmon satellites or exciton-exciton coupling. Our approach starts from the GW + Bethe-Salpeter approximation to many-body perturbation theory which is based on a quasiparticle picture, and it adds coupling of excitations through a consistent inclusion of dynamically screened interactions. This requires one to consider scattering contributions that are usually neglected. The result is formulated in a way that highlights essential physics, that can be implemented as a postprocessing tool in first-principles codes, and that suggests which kind of materials and measurements should exhibit strong effects. This is illustrated using a model.

The quantum-mechanical nature of electrons is responsible for many striking properties of materials, such as the fact that light absorption can happen only at certain wavelengths. On top of this, the Coulomb interaction between all the particles can drastically influence these observations. In particular, photoemission or optical absorption spectra are in general very different from what one would expect from noninteracting electrons, exhibiting shifts and broadening of peaks, as well as additional structures called satellites [1]. Though often difficult to interpret, they are tightly linked to technologically important phenomena such as exciton relaxation and decoherence, multiple exciton generation [2][3][4], or singlet exciton fission observed in molecular crystals [5]. The underlying mechanism is in all cases a coupling of excitations.
In the framework of many-body perturbation theory [6,7] this coupling can be expressed through the dynamically screened Coulomb interaction W , where the bare interaction is screened by the charge response. To lowest order in W electron addition and removal spectra are obtained in the widely used GW approximation for the electron self-energy [8]. This reflects a picture of a quasiparticle coupled to neutral charge excitations, which create satellites because W (ω) is frequency dependent. The spectrum of these excitations, and hence W , can in turn be obtained by solving the two-particle Bethe-Salpeter equation (BSE) for the irreducible polarizability P, where the dressed electrons and holes interact to lowest order through the dynamically screened W itself [7,9]. In this way neutral excitations are coupled to each other, which can lead to satellites also in the spectra of P or W , and to multiexcitation effects such as those mentioned above.
However, most state-of-the-art BSE calculations neglect the frequency dependence of W [1,9]. They can therefore not access important coupling effects such as the damping of quasiparticles and double plasmon excitations observed in inelastic x-ray spectroscopy of aluminum and sodium at large momentum transfer [10][11][12][13], or multiple excitations in openshell molecules [14], in closed-shell systems such as polyene [15], or in strongly correlated materials such as NiO [16,17]. For real materials, only few calculations beyond the static BSE have been performed [7,[18][19][20][21][22][23], mostly for quasiparticle excitation energy, lifetime, and weight renormalization, while satellites were not addressed. The most studied situation is core-level absorption, where quasiboson models led to important insight [24][25][26][27][28]. Core-level absorption satellites were also described from first principles but making drastic approximations, in particular concerning static excitonic effects [29]. Double excitations have been studied in model systems (see, e.g., [30][31][32][33]) including double plasmon satellites in the homogeneous electron gas or simple metals [10,11,34,35]. A generally applicable approach is instead, to the best of our knowledge, still missing.
The fully dynamical BSE is complicated [1], and its solution might not be worth the effort since the BSE is similar to the Dyson equation for the one-body Green's function G in the GW approximation, which often fails for satellites [1,[36][37][38][39][40][41]. Instead, a cumulant approach [42,43] for P, reflecting a picture of coupled bosons, is more promising. It is the exact solution for a two-level limiting case [44], and it is additionally motivated by the success of an increasing number of ab initio calculations for the one-body G using a cumulant of first order in W [37,41,[44][45][46] while higher order contributions appear to be negligible [47][48][49]. The aim of the present work is to derive a cumulant approximation for the polarizability with similar computational requirements as a static BSE calculation and applicable to core and valence excitations of a broad range of materials, in order to open the way for the understanding and prediction of coupling effects in realistic materials.
We start with the known perturbation expansion for P in terms of W and the Hartree Green's function G H [50]. Usually such an expansion is partially resummed through a Dyson equation, which in the present case would be the BSE. Here we want to achieve a better representation of the perturbation series. The challenge is to arrange and approximate the perturbation expansion in a way that suggests a cumulant representation. We start by considering the contributions up to second order in W ; all typical second-order diagrams are shown in Fig. 1. In state-of-the-art Bethe-Salpeter calcula- , an instantaneous screened Coulomb interaction. In this case the five diagrams of type 1(f)-1(j) vanish in the Tamm-Dancoff approximation (TDA) [1,51,52], where elementary scattering processes between empty and occupied states are neglected. The remaining diagrams taken to infinite order yield the static BSE (SBSE) P(1423) =P 0 (1423) +P 0 (1221)W 0 (12)P(1423). (1) Here, 1 ≡ (x 1 t 1 ) ≡ (r 1 σ 1 t 1 ) stands for space, spin, and time and f (n) ≡ dn f (n). Moreover,P 0 (1423) ≡Ḡ (13) GW self-energy where W → W 0 [8,53]. In principle the full BSE for P(1423) with four space, spin, and time arguments has to be solved to derive P(1313), which yields electron-hole excitation spectra. However, the instantaneous W 0 allows one to solve Eq. (1) directly in two times, or one frequency.
To overcome the static approximation, we consider corrections to first order in W ≡ W − W 0 , stemming from first-order diagrams in W and beyond, which are not zero in TDA. The second-order diagrams 1(a)-1(e) yield corrections that are first order in W 0 and first order in W , but now also 1(g) contributes as well as 1(h) when W 0 is an electron-hole interaction (vertical interaction line) and W dresses the Green's functions (horizontal). Instead, when in 1(h) W 0 is the dressing interaction and W contributes to the electron-hole interaction, the diagram vanishes in the TDA. This analysis is important, because the nonvanishing diagrams and all their higher orders in W 0 simply transform G H intō G and add instantaneous electron-hole interaction (ladder) lines. To first order in W one obtains the two prototypical diagrams in Fig. 2, in which arrows represent now dressed Green's functionsḠ, dashed lines stand for W 0 , wiggly lines for W , and the red features should not be considered at this stage. Since our choices are based on the TDA, the approach is valid for semiconductors and insulators at low temperature. For metals or strongly doped semiconductors one may have to take into account diagrams 1(i)-1(j), which express changes in the dynamical screening of the system due to an excitation. This case has been addressed in Ref. [54].
Our next goal is to express the diagrams in Fig. 2 in terms of W and the two-times solution of the SBSE, P(x 1 t 1 x 4 t 3 x 2 t 1 x 3 t 3 ). It is impossible to detectP in the black diagrams of Fig. 2, but using the relation dx 2 Ḡ (x 2 t 1 2)Ḡ(2 x 2 t 1 ) = ±iḠ(2 2) (t 2 t 1 t 2 ), (2) (with the + and − signs for hole and electron propagators, respectively) which holds when t 1 lies between t 2 and t 2 as indicated by (t 2 t 1 t 2 ), allows us to insert additional spacespin points that are integrated over (red dots in Fig. 2). Their time coordinate can be chosen equal to an already existing time integration point, as indicated by the red ovals. Now, taking t 2 = t 1 and t 4 = t 3 , the figure clearly exhibits the lesser part (i.e., t 1 < t 3 ) of the two-times polarizability P < in terms ofP and W (equivalently for P > , when t 1 > t 3 ). In a pair basis labeled by λ this reads where t 13 stands for (t 1 , t 3 ), repeated indices are summed over, and λλ (t 13 ) given by is an effective exciton self-energy, becauseP λλ (t 13 ) represents the propagator of independent excitons. It is the analog of the GW self-energy for one-particle excitations with G and W replaced byP and W, respectively. The effective excitonexciton interaction W is given by matrix elements of W that can be read from Fig. 2 and that are detailed in [55]. W = W pp + W eh consists of an effective electron-electron or hole-hole interaction labeled pp from the first diagram in Fig. 2, and an electron-hole (eh) interaction from the second diagram. This reflects the fact that excitons are composite particles [56], and their effective interaction results from the interaction between its constituents, i.e., all electrons and holes. The two terms have opposite signs, leading to partial cancellation of dynamical effects as suggested in Ref. [20] and discussed below. This compact result cannot be obtained with the usual derivation of the BSE from the GW approximation, since it necessitates all the nonvanishing diagrams in Fig. 1.
In the excitonic basis (solution of the SBSE), which mixes transitions at different k points in the Brillouin zone,P is diagonal, withP < λλ (t 13 ) = e −iE λ (t 3 −t 1 ) . Here E λ is the excitation energy, and λ = (λ, q) labels specifically the exciton band index λ and the exciton wave vector q. Moreover, we suppose that the exciton self-energy in (3) is diagonal in this basis. Comparison of the resulting diagonal of Eq. (3) to the firstorder expansion of a cumulant ansatz In the limit of instantaneous interaction, this yields the standard SBSE result, such that static excitonic effects are fully taken into account. The cumulant expression can be understood in terms of a multilevel system coupled to bosonic modes. Our result fulfills the TDA limit of a two-level system coupled to one boson mode that was derived in Ref. [44]. The bosonic modes also contain excitons when W 0 and W are calculated in the SBSE. Finally, for a core hole and when excitonic effects in the SBSE are negligible, Eq. (6) reduces to the particle-core-hole cumulant derived in Ref. [29]. The polarizability resulting from (5) and (6) can be written in terms of a renormalization factor e −R λ , a correction E λ with respect to the SBSE exciton energy, and a termC λ responsible for satellite structures: where , analogous to the one-body cumulant G [37]. For a model system consisting of two free electron bands coupled with a nondispersive boson of frequency ω 0 through a coupling constant g, and surrounded by a homogeneous dielectric with dielectric constant , the solutions of the SBSE are Wannier excitons [57] with energy E λ = QP − μ , where QP is the quasiparticle (QP) band gap, μ is the reduced mass of the electron-hole pair, n λ is the Rydberg quantum number of state λ, M is the exciton mass, and q the modulus of its wave vector [55]. Figure 3(a) shows the SBSE and cumulant spectral functions of P at q → 0, which corresponds to an absorption spectrum, for the lowest excited state (Wannier exciton in the 1S configuration; we omit the subscript λ in the following). The dynamical effects create a slight redshift E of the exciton QP peak, a series of satellites at energies E + E + nω 0 with integer n, and a weight transfer from the QP peak to the satellites, due to e −R . Figure 3(b) highlights pp, eh, and (for the second satellite) mixed contributions. In the first satellite pp and eh have opposite signs, but the negative eh contributions are always smaller. Since the resulting structures are located at the same frequency, the spectral function remains positive. The strength of the second satellite goes as W 2 according to an expansion of Eq. (5). Therefore only the mixed contribution stemming from the product pp-eh is negative. The weight transfer from the QP peaks to the satellites quantified by e −R depends for fixed coupling strength g on and on the ratio between the exciton bandwidth and the boson frequency, d = ω 0 , where q 0 is the radius of the first Brillouin zone, assumed to be spherical. Figure 4 shows the renormalization factor e −R , as a function of the exciton Bohr radius and of the ratio d between the exciton band-width and the boson frequency. The weight transfer is important (e −R 1) for nondispersive excitons (small d), and much weaker when d 1. Indeed, the exciton bandwidth is related to the rate of the exciton hopping processes while the boson frequency ω 0 expresses the timescale of the charge fluctuations associated to the excitations (excitons, free electron-hole pairs, or plasmons) of the system. When d 1 the hopping processes are fast compared to the charge fluctuations and the exciton behaves like a particle propagating in a static medium, so W is negligible. For fixed d, the evolution of e −R with π a/r 0 shows that a strongly localized exciton is much less influenced by dynamical effects than a delocalized one. This can be attributed to the fact that R ≡ R pp + R eh is a sum of pp (R pp ) and eh (R eh ) contributions, with R pp positive and R eh negative and vanishing for a r 0 → ∞. Therefore R pp and R eh cancel to a large extent for small a, which is intuitive: for two electron-hole pairs A and B placed with the same relative coordinates in two adjacent unit cells, the electron (hole) of A and the electron (hole) of B are found at a distance of the order of r 0 , independently of the exciton radius. The distance between the electron (hole) of A and the hole (electron) of B, instead, is of the order of r 2 0 + a 2 , which is always larger. This explains why |W eh | < |W pp | for each a. For a r 0 , the eh distance is proportional to a, independently of r 0 , such that |W eh | becomes negligible; this trend will be enhanced by the fact that screening is stronger at larger distances. When a r 0 the typical length scale is r 0 , the same for W pp and W eh , and the cancellation becomes exact. Of course, in real systems also the different character of the valence and conduction states makes the cancellation imperfect, but this will not change the trends. Our results suggest therefore that dynamical effects are very sensitive to the degree of localization of the exciton, and that for a strongly localized exciton the pp and the eh dynamical effects cancel, which explains Fig. 4. Screening plays a key role here: on one side, there is an overall scaling of the R's with the inverse dielectric constant, meaning that weak screening makes the single terms larger, which may overall increase dynamical effects. On the other hand, the exciton Bohr radius is smaller when screening is weaker, which favors cancellations and therefore a decrease. In practice, satellites in core and valence absorption or in loss spectra are much less observed than in photoemission spectra, pointing to the fact that the cancellation effects play a predominant role. Since they depend on the exciton localization, they can only be captured correctly when static excitonic effects are taken into account, as is the case in our formulation.
In conclusion, we have derived a cumulant formulation for excitation spectra that contains excitonic effects and the coupling between excitons or other neutral excitations. It uses as input the results of a standard GW + BSE calculation where the frequency dependence of the screened Coulomb interaction is neglected, and it adds dynamical effects in a simple, transparent, and easy to implement way. In order to obtain a consistent expression, diagrams had to be taken into account that can be neglected in the absence of dynamical effects. In some important limiting cases, our formulation reduces to approximations found in the literature. We have applied the approach to a simple model, which allows us to highlight cancellations that justify the use of the static BSE in certain parameter ranges, and to estimate other regions in parameter space where materials should instead exhibit sizable dynamical coupling effects. This information should be precious in the search for materials where these effects could be used for applications.