# The valence electron photoemission spectrum of semiconductors: ab initio description of multiple satellites

###### Abstract

The experimental valence band photoemission spectrum of semiconductors exhibits multiple satellites that cannot be described by the GW approximation for the self-energy in the framework of many-body perturbation theory. Taking silicon as a prototypical example, we compare experimental high energy photoemission spectra with GW calculations and analyze the origin of the GW failure. We then propose an approximation to the functional differential equation that determines the exact one-body Green’s function, whose solution has an exponential form. This yields a calculated spectrum, including cross sections, secondary electrons, and an estimate for extrinsic and interference effects, in excellent agreement with experiment. Our result can be recast as a dynamical vertex correction beyond GW, giving hints for further developments.

Photoemission is a prominent tool to access information about electronic structure and excitations in materials. Modern synchrotron sources can provide detailed insight, thanks to their high intensity and broad photon energy range. But the interpretation of the experimental data is far from obvious, and theory is an essential complementary tool. However, ab initio calculations typically focus on bulk bandstructure rep.prog.phys.61.237 (); ssp.54.1 (); thus surface effects are ignored, and satellites are not included. The latter are a pure many-body effect due to coupling to excitations of the material. Such many-body effects are contained in approaches developed for correlated materials georges96 (); kotliar06 () however, these are usually based on models with short-range interactions, whereas satellites such as plasmons involve long-range effects. Plasmon satellites have been extensively studied in core-level experiments [e.g.\@:][]PhysRevB.76.085422. There they can be described by a theoretical model where a single dispersionless fermion couples to bosons. The resulting exact Green’s function has an exponential form given by the so-called cumulant expansion (CE). A Taylor expansion of the exponential leads to a well defined quasi-particle (QP) peak followed by a decaying series of plasmon satellites at energy differences given by the plasmon energy, consistent with experimental observations HEDIN80 (); Hedin99 (); NOZIERES69 (); LANGRETH70 (); Bechstedt2003 (). In the valence region, plasmon satellites are much less studied, though ab initio approaches can provide a good starting point. At high photoelectron energies the photoemission spectrum is approximately proportional to the intrinsic spectral function , where is the one-particle Green’s function. The latter is typically calculated using the widely used GW approximation (GWA) Hedin65 (); Hedin69 (); Hedin99 (). In principle, the GWA contains correlations effects beyond the quasiparticle approximation. However, these additional features are rarely calculated due to computational complexity and, more importantly, the serious discrepancies between GWA and experiment (see e.g. PhysRevB.46.13051 (); Aryasetiawan96 (); ISI:000179080800133 (); Kheifets03 ()). The CE has also been used for the homogeneous electron gas Holm97 () and simple metals ISI:000179080800133 (); Aryasetiawan96 (), yielding an improved description of satellites over GW. Silicon Kheifets03 () and graphite ISI:000166608600008 () were also studied, but no plasmon satellite series were observed. However, these results are not conclusive due to difficulties of interpreting the experimental data. This leaves a series of important questions: (i) do materials generally exhibit intrinsic satellites in the valence band region following a cumulant like distribution, or are the extrinsic plasmon peaks, due to losses incurred by the escaping photoelectron, dominant? (ii) if such series are seen, how bad are ab initio GW calculations, what is the reason for their failure, and (iii) how can they be improved? Answering these questions would be a crucial step towards a better understanding of correlation effects in electronic excitations and a predictive ab initio approach to photoemission.

In this work we focus on plasmon satellites using silicon as a prototypical example. We have obtained valence band photoemission data at high photon energy (XPS) that constitute a reliable and well resolved benchmark. Analysis of the data allows us to elucidate the failure of GW in describing the satellites. Then, starting from the fundamental equations of many-body perturbation theory (MBPT), we show how the failure can be overcome by using a decoupling approximation that leads to an exponential representation of the one-particle Green’s function. Together with an estimate for extrinsic and interference effects, we obtain results for the quasiparticle peaks and satellites in excellent agreement with experiment. Our theoretical results can be expressed in terms of a dynamical vertex correction, a powerful basis for further modelling.

Angular resolved valence photoemission (ARPES) measurements were performed at the UHV photoemission experimental station of the TEMPO beamline exp1 () at the SOLEIL synchrotron radiation source. Linearly polarized photons from the Apple II type Insertion Device (HU44) were selected in energy using a high resolution plane grating monochromator with a resolving power . The end-station chamber (base pressure mbar) is equipped with a modified SCIENTA-200 electron analyzer with a delay-line 2D detector which optimizes the detection linearity and signal/ background ratio ISI:000287530900019 (). The overall energy resolution was better than 200 meV. The photon beam impinges on the sample at an angle of 43, and photoelectrons were detected around the sample surface normal with an angular acceptance of . A -type ( atoms/cm) Si(001) wafer was cleaned from the native oxide by flash annealing at 1100 C after prolonged degassing at 600 C in ultra-high vacuum. The silicon surface was annealed at 300 C to prevent surface etching, and hydrogenated in a partial pressure of activated hydrogen about mbar for 20 min. The ARPES was measured along the direction. At 800 eV kinetic energy the Si Brillouin zone is observed with an emission angle slightly smaller than 5. The measured photoemission map was integrated over the spectral intensity originated by two Brillouin zones. The Fermi level was obtained by measuring a clean Au(111) surface.

The experimental data (crosses) are summarized in Fig. 1. One can distinguish the quasiparticle peaks between the Fermi level at zero and the bottom valence at -12 eV, followed by two prominent satellite structures, each at a mutual distance of about 17 eV, as well as a more weakly visible third satellite between -52 and -60 eV. These structures are obviously related to the 17 eV silicon bulk plasmon ISI:A1978FY99000003 ().

The exact one-electron Green’s function is described by an equation of motion with the form of a functional differential equation kadanoffbaym (),

(1) |

Here is the non-interacting Green’s function, is a fictitious external perturbation that is set to zero at the end of the derivation, is the bare Coulomb interaction, and all quantities are understood to be matrices in space, spin, and time. The Hartree potential gives rise to screening to all orders. Linearizing with respect to yields 2011arXiv1103.1630L ()

(2) |

where is equal to screened by the inverse dielectric function, is the screened Coulomb interaction, and is the Green’s function containing the Hartree potential at vanishing ; only time arguments are displayed explicitly and repeated indices are integrated. This linearization preserves the main effects of and hence of plasmons. With the additional approximation one obtains the Dyson equation in the GWA for the self-energy . However this approximation can be problematic. For the following analysis we use the standard approach, where is taken from an LDA calculation and is the screened interaction in the Random Phase Approximation.

Fig. 2 shows the spectral function Fleszar97 () of Si at the point, for top valence (solid line) and bottom valence (dashed), respectively. The top valence shows a sharp quasiparticle peak followed by a broad, weak satellite structure at about -21 eV. This peak stems from the prominent peak in (full circles) at about eV, itself due to the plasmon peak in . It is a typical plasmon satellite, though (cf. Hedin99 ()), the QP-satellite spacing is slightly overestimated because the term (full squares) in the denominator of the expression for is not constant. However the GWA has a more severe problem: for the bottom valence, the satellite structure at about eV is much too far from the QP peak at about eV, and much too sharp. This satellite does not correspond to a plasmon peak in (empty circles), but to a zero in (empty squares) in the denominator of , as for a QP peak. It has been interpreted in the HEG as a plasmaron, a coupled hole-plasmon mode Lundqvist1967a (); *ISI:A1968A818100003, but as noted below it is an artifact of the GWA ISI:A1972N855800007 (); ISI:A1973P093700012 (). Fig. 1 compares the total spectral function (dashed red line) summed over all valence bands and -points, with our XPS data. The effects of cross-sections are included by projecting on angular momenta in atomic spheres using the atomic data of Ref. Yeh85 (). The secondary electrons background at energy was determined by integrating the calculated intrinsic spectral intensity between and the Fermi level, similar to PhysRevB.5.4709 (). A constant scaling factor was set such that the measured photoemission intensity at the highest binding energy (60 eV), where primary electrons intensity is absent, is reproduced. As expected, the dominant QP spectrum is well described by , but the satellite is dominated by the plasmaron around eV, in complete disagreement with experiment. The experimental plasmon satellite at about eV appears only as a weak shoulder in the . Thus the plasmaron peak is responsible for the GWA failure ISI:A1972N855800007 (); ISI:A1973P093700012 () in silicon.

To go beyond the GW self-energy requires vertex corrections. However, adiabatic vertex corrections (see e.g. PhysRevB.49.8024 ()) only lead to renormalization of energies and do not create new structures. Thus alternatively, we concentrate here on dynamical effects, and we choose to approximate directly Eq. (2), without passing through a self-energy. We decouple Eq. (2) approximately by supposing that and are diagonal in the same single particle basis. Eq. (2) is then applicable separately for every single matrix element of and each state couples independently to the neutral excitations of the system through [InasimilarspirittheequationofmotionfortheGreen'sfunctionisdecoupledin:][]Almbladh83. The latter can now be understood as the screened intra-orbital Coulomb matrix element for the chosen state. Such a decoupling approximation can be optimized ISI:A1972N855800007 (); ISI:A1973P093700012 () by adding and subtracting a self-energy correction, hence by using a QP Green’s function obtained from a good QP self-energy instead of . Since the GWA is currently the state-of-the art for QP properties, we suppose that for every decoupled state , is determined from , where is the (complex) GW quasiparticle energy and . Now Eq. (2) can be solved exactly for each state. Briefly the main steps are: (i) solve the non-interacting () version of (2), which leads to an explicit solution ; (ii) iterate the result starting from . Here compensates for the self-energy insertion used for the optimized decoupling; (iii) use the exact relation to derive

(3) |

The equilibrium solution is obtained setting .

In silicon, where the peaks in the loss function are well defined, it is justified to use a single plasmon pole model with plasmon energy and intrinsic strengths for each matrix element of . Besides , the total exponent becomes then with obtained from the corresponding GW results. We find that varies around 0.3. Taylor expansion of the exponential leads then to the spectral function

(4) |

where and . Eq. (4) is similar to the plasmon pole version of the CE (cf. Ref. Aryasetiawan96 ()). However here the exponential solution arises from a straightforward approximation to the fundamental differential equation (1): the linearization of the Hartree potential reveals the boson of the model (i.e., the plasmon via peaks in ), and the diagonal approximation of gives rise to each isolated fermion. Our results are summarized in Fig. 1. The dot-dashed line gives the result of this procedure together with the cross sections and the secondary electron background. The shapes of the QP peaks change little with respect to GW, but now the full series of satellites is present. The internal structure of the satellites which originate from the multiple valence bands, is also reproduced. This validates the decoupling approximation in the dense valence band region where, contrary to the case of an isolated core level, its success is a priori far from obvious. However, the intensity of the observed satellites is significantly underestimated. This discrepancy is similar to that found for the CE in simple metals, where extrinsic losses were suggested as a likely cause Aryasetiawan96 (). These might also be reduced by interference effects ISI:A1997WK15700033 (). To check this possibility we estimated the contributions from both effects to the satellite strengths using Eq. (32) and (36) of Ref. HedinMichielsInglesfield98 (). This approach uses a plasmon pole model, Inglesfield fluctuation potentials, and an average over hole positions that takes account of the photoelectron mean free path HedinMichielsInglesfield98 (). We observe that the averaged total satellite line shape in this model is similar to that for the intrinsic part, with a width eV due to plasmon dispersion. Thus we can approximate the extrinsic and interference effects by renormalizing the intrinsic satellite intensity, i.e., by the replacement in Eq. (4). These quantities are evaluated with eV and Å at 800 eV for Si, yielding = 0.63 and . This also modifies the strength of the QP peaks, but preserves overall normalization. The broadening of the satellites must also be increased, . The total spectrum thus obtained (black line) is in unprecedented agreement with experiment. We stress that this result contains no fit parameters besides the two scaling factors (for spectrum and background) due to the arbitrary units of the experiment.

The success of our present approach stresses the need to go beyond the GWA. The exponential representation of implicitly corresponds to a vertex correction to the self-energy. Since our derivation yields as a function of the screened potential (3), this functional derivative can be performed explicitly, using . From Eq. (3), a straightforward derivative of contains a series of satellite contributions. The two inverse Green’s functions lead to a significant complication, because they contain the inverse of this series. This clearly illustrates the difficulty of modelling in order to treat dynamical effects. It suggests rather to concentrate on modelling , where the various contributions are simply summed, and hence to search for a self-energy of the form instead of . In conclusion, on the basis of our experimental XPS data we have analyzed the failure of GW to reproduce plasmon satellites and linked this failure to the appearance of a artificial plasmaron peak. On the other hand, GW results are fair when the imaginary part of , hence the intensity of the corresponding plasmon, is small enough so that no sharp plasmaron is created. Thus surprisingly, one might expect GW to work better in describing satellites stemming from local plasmon or interband excitations close to the Fermi level in “strongly correlated” materials than for the strong plasmon structures in conventional semiconductors. Starting from the fundamental equations of MBPT we have derived an exponential solution to the one-particle Green’s function, analogous to that from the CE, that overcomes the drawbacks of the GWA. Comparison to new photoemission data shows that this yields a very good description of the spectral function of bulk silicon, including the satellites series. By calculating the secondary electron background, cross section corrections as well as a correction for extrinsic and interference effects, we achieve an agreement between theory and experiment that can be considered as a benchmark. Our derivation also suggests how the results can be improved in cases where the presently used approximations are inadequate. Finally, by accessing an expression for the vertex function, our approach yields precious hints for directions to take in modeling dynamical effects beyond the GWA.

###### Acknowledgements.

We acknowledge support by ANR (NT09-610745), St Gobain R&D (091986), and Triangle de la Physique (2007-71). JJR and JJK are also supported in part by DOE BES Grant DE-FG03-97ER45623 and DOE CMCSN. Computer time was granted by IDRIS (544). We thank F. Bechstedt for fruitful discussions.## References

- (1) F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61, 237 (1998)
- (2) W. Aulbur, L. Jonsson, and J. Wilkins, Solid State Physics 54, 1 (2000)
- (3) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996)
- (4) G. Kotliar et al., Rev. Mod. Phys. 78, 865 (2006)
- (5) F. Offi et al., Phys. Rev. B 76, 085422 (2007)
- (6) L. Hedin, Phys. Scripta 21, 477 (1980)
- (7) L. Hedin, J. Phys.: Cond. Mat. 11, R489 (1999)
- (8) P. Nozières and C. De Dominicis, Phys. Rev. 178, 1097 (1969)
- (9) D. Langreth, Phys. Rev. B 1, 471 (1970)
- (10) F. Bechstedt, “Principles of surface physics,” (Springer, 2003) Chap. 5, pp. 199–201
- (11) L. Hedin, Phys. Rev. 139, A796 (1965)
- (12) L. Hedin and S. Lundqvist, Solid State Physics, Vol. 23 (Academic Press, New York, 1969)
- (13) F. Aryasetiawan, Phys. Rev. B 46, 13051 (1992)
- (14) F. Aryasetiawan, L. Hedin, and K. Karlsson, Phys. Rev. Lett. 77, 2268 (1996)
- (15) M. Vos et al., Phys. Rev. B 66, 155414 (2002)
- (16) A. S. Kheifets et al., Phys. Rev. B 68 (2003)
- (17) B. Holm and F. Aryasetiawan, Phys. Rev. B 56, 12825 (1997)
- (18) M. Vos, A. S. Kheifets, E. Weigold, and F. Aryasetiawan, Phys. Rev. B 63, 033108 (2001)
- (19) F. Polack et al., AIP conf. proc. 1234, 185 (2011)
- (20) N. Bergeard et al., J. Synch. Rad. 18, 245 (2011)
- (21) J. Stiebling, Zeitschrift für Physik B-Condensed Matter 31, 355 (1978)
- (22) L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics (W.A. Benjamin Inc., New York, 1964)
- (23) G. Lani, P. Romaniello, and L. Reining, ArXiv e-prints(2011), arXiv:1103.1630 [cond-mat.str-el]
- (24) A. Fleszar and W. Hanke, Phys. Rev. B 56, 10228 (1997)
- (25) B. Lundqvist, Zeitschrift für Phys. B Cond. Mat. 6, 193 (1967)
- (26) B. Lundqvist, Phys. Kondens. Materie 7, 117 (1968)
- (27) C. Blomberg and B. Bergerse, Can J. Phys. 50, 2286 (1972)
- (28) B. Bergerse, F. Kus, and C. Blomberg, Can J. Phys. 51, 102 (1973)
- (29) J. Yeh and I. Lindau, Atom. Data Nucl. Data Tables 32, 1 (1985)
- (30) D. A. Shirley, Phys. Rev. B 5, 4709 (1972)
- (31) R. Del Sole, L. Reining, and R. W. Godby, Phys. Rev. B 49, 8024 (1994)
- (32) C.-O. Almbladh and L. Hedin, “Handbook on synchrotron radiation,” (North-Holland, Amsterdam, 1983) Chap. 8
- (33) F. Bechstedt, K. Tenelsen, B. Adolph, and R. DelSole, Phys. Rev. Lett. 78, 1528 (1997)
- (34) L. Hedin, J. Michiels, and J. Inglesfield, Phys. Rev. B 58, 15565 (1998)