Modeling core-hole screening in jellium clusters using density functional theory

The screening of a 2p core-hole in Na clusters is investigated using density functional theory applied to an extended jellium model with an all-electron atom in the center. The study is related to recent experiments at the free electron laser at DESY in which photoelectron spectra from mass-selected, core-shell ionized metal clusters have been recorded. Relaxed and unrelaxed binding energies as well as Kohn-Sham orbital energies are calculated in Perdew-Zunger self-interaction-corrected exchange-only local spin-density approximation for valence and 2p core electrons in Na clusters up to 58 atoms. The relaxed binding energies follow approximately the metal-sphere behavior. The same behavior is seen in the experiment for sufficiently big clusters, indicating perfect screening and that the relaxation energy due to screening goes to the photoelectron. Instead, calculating the kinetic energy of the photoelectrons using unrelaxed binding energies or Kohn-Sham orbital energies yields wrong results for core-shell electrons. The screening {\em dynamics} are investigated using time-dependent density functional theory. It is shown that screening occurs on two time scales, a core-shell-dependent inner-atomic and an inter-atomic valence electron time scale. In the case of Na 2p ionization the remaining electrons in the 2p shell screen within tens of attoseconds while the screening due to cluster valence electrons occurs within several hundreds of attoseconds. The screening time-scales may be compared to the photon energy and cluster size-dependent escape times of the photoelectron in order to estimate whether the photoelectron is capable of picking up the relaxation energy or whether the residual system is left in an excited state.


Introduction
Core-holes are created when matter is irradiated with energetic photons. The net binding energy is obtained from the kinetic energy of the photoelectron upon subtracting the photon energy. As the target size increases from the isolated atom via molecules and clusters to, ultimately, bulk matter, the shift of the net binding energy contains information about the screening capability of the other electrons in the system and allows us to investigate, e.g., metal-to-nonmetal transitions. In fact, x-ray photoelectron spectroscopy (XPS) is well established and extensive work has been devoted to the measurement of core-level binding energy shifts in solids and at surfaces [1]. However, the systematic study of core-level shifts as a function of the target size became possible only recently because of the necessity of having size-selected clusters and tunable, powerful x-ray sources such as free-electron lasers available.
From the theoretical point of view, the calculation of core-level binding energies in manyelectron systems is nontrivial, both conceptually and computationally [1,2]. The study of atoms, clusters or bulk, in which suddenly a core-electron is removed by single-photon ionization, has a long history, starting probably from Slater's transition state theory [3,4] and various kinds of self-consistent field methods ( SCFs) [5,6]. SCF approaches are based on energy balances between electronic configurations, and it is not always clear whether relaxed or unrelaxed configurations should be taken. Here, 'relaxed' means that the other electrons are allowed to assume the new, energetically most favorable configuration-albeit without filling the core-hole. 'Unrelaxed' means that the other electrons are kept frozen. It is intuitively clear that the unrelaxed energy difference between the final and the initial configuration should be considered if the photoelectron escapes too rapidly to notice the relaxation dynamics of the other electrons [7][8][9].
Methodologically, self-consistent field methods such as Hartree-Fock (HF) or density functional theory (DFT) (see, e.g., [10]) are usually employed for the calculation of energy differences between the final and the initial configuration in SCF approaches. As both methods introduce orbitals and orbital energies, the interesting question arises of how these orbital energies are related to the (relaxed or unrelaxed) binding energies [6,11,12]. Clearly, the energy difference depends on the initial and the final state, while the orbital energies do not, unless fractional occupation numbers, as in Slater's transition state theory, are introduced. 3 We use an extended Na-cluster jellium model with an all-electron Na atom in the center in order to have core-levels in the system at all. This model may be viewed as a hybrid of the pure cluster jellium model [13] and bulk-models in which atoms are immersed in a homogeneous electron gas [5,14]. We employ DFT to calculate the relaxed and unrelaxed energy differences between the initial ground state and the final configuration with either a valence electron or a 2pcore electron removed. The results as a function of the cluster size (up to 58 atoms) are compared with the prediction of the metal-sphere model. Time-dependent DFT (TDDFT) (see [15] for a state-of-the-art account) is used to investigate the screening dynamics after the sudden removal of a 2p electron. We apply our model in this first study to Na clusters (instead of Pb − clusters in the experiments [16,17]) in order to keep the number of inner electrons in the embedded atom and thus the numerical effort manageable.
This paper is organized as follows. The basic methodology and theoretical background are introduced in section 2. The DFT results on the size dependence of binding energies and the TDDFT results on the screening dynamics are presented in section 3. We conclude in section 4.

Basic theory and models
We use DFT in exchange-only local spin-density approximation (xLSD) with the Perdew-Zunger (PZ) self-interaction correction (SIC) [11]. The electronic Kohn-Sham (KS) equation reads Here, iσ is the energy of KS orbital ϕ iσ , i = 1, 2, . . . , N σ with N σ , being the number of electrons with spin-projection σ =↑, ↓. T is the kinetic energy operator and V (r ) is the spherically symmetrical external potential introduced below. V H and V xcσ are the Hartree potential and the exchange-correlation (xc) potential in xLSD approximation, respectively (see, e.g., [10]). Both are functionals of the (spin) density. The term in curly brackets is PZ-SIC; that is, in the KS equation for orbital ϕ iσ , the Hartree-xc potential V Hiσ + V xciσ evaluated with the corresponding spin density n iσ = |ϕ iσ | 2 is subtracted from the full Hartree-xc potential It is known that the PZ-SIC applied to xLSD significantly improves the results as compared to uncorrected xLSD [11,19]. The central-field approximation is applied to V H + V xcσ − {V Hiσ + V xciσ } in (1) so that the single-particle orbital angular momentum quantum numbers and magnetic quantum numbers m remain good quantum numbers. Although the PZ-SIC in general improves the uncorrected KS results quantitatively, there have been several critical issues discussed in the literature [18,19], one of them being the orbital-dependent effective Hamiltonian in (1), leading to non-orthogonal KS orbitals. In our calculations we therefore force the KS orbitals to be orthogonal. An extended version of the Qprop code [20] (http://www.qprop.de) is used to solve (1) and-for the dynamics in section 3.2-its timedependent version (in which iσ is replaced by i∂ t and the respective adiabatically timedependent potentials are used [15]).

Cluster jellium model with the central all-electron atom
The external potential V (r ) in (1) of this work is The first term −Z /r is the pure Coulomb potential of the central atom for which all electrons are taken into account. This atom provides the ionizable core-shell. The other A − 1 atoms in the cluster are treated in jellium approximation. The jellium potential is that of a homogeneously charged sphere of radius mimicking the net ionic background of the cluster, as 'seen' by the cluster valence electrons. A Wigner-Seitz radius for Na of r s = 4 was used. For A = 1 the potential (3) clearly reduces to the single-atom case. Figure 1 shows the results for electron density, the PZ-KS potentials and the KS levels for the neutral closed-shell Na 58 cluster. Because all electrons in the central atom are considered 5 there is a peak in the electron density in the center. What would be the 3s level in the isolated Na atom is now the 1s jellium cluster level. The higher the jellium level the closer the respective KS eigenvalue is to the one in the corresponding pure valence electron jellium model. Note that the PZ-SIC yields the correct 1/r -behavior of the KS potentials, unlike pure xLSD.

Determination of the core-shell binding energies
The experimental observable is the kinetic energy kin of photoelectrons that are emitted upon irradiation of clusters with x-rays of photon energyhω. Energy conservation requires that so that the binding energy can be determined as E 0 < 0 is the initial energy of the target. In what follows we assume that the target is initially in its ground state. However, it is far from obvious what E f is. As long as Auger decay occurs on time scales longer than the photoelectron escape times, the final target configuration of interest to us is a cluster in which one atom has a core-hole 1 . We assume that the photon energy is such that most likely one of the 2p electrons is removed. If the escape time is longer than the screening time the photoelectron may pick up the relaxation energy and the net binding energy is On the other hand, if the photoelectron escapes before the system is relaxed, the binding energy reads According to Koopmans' theorem the unrelaxed energy difference is equal to the Hartree-Fock orbital energy, − HF = E HF unrlxd − E HF 0 > 0. The fact that in the ultra-high photon energy limit the centroid of the photoelectron spectrum is located at kin = HF +hω is known as the Manne-Åberg theorem [7] or the Lundqvist sum rule [8]. However, Koopmans' theorem does not hold in DFT. Instead, the 'non-Koopmans energy' tends to cancel a part of the relaxation energy (7) [11] so that KS orbital energies | | = − are shifted towards E rlxd . Whether − is above or below E rlxd depends on the xc energy functional used. Introducing fractional occupations 0 f 1, Janak's theorem establishes a connection between the relaxed binding energy and the relaxed KS orbital energies ( f ) (see, e.g., [11]).
Using the extended cluster jellium model introduced in section 2.1, we determined the binding energy of a core-shell electron iσ by switching the fractional occupation number of the respective KS orbital f iσ from 1 to 0. Without relaxing the final configuration with the core-hole, we obtain 6 and with relaxation The relaxation energy in general, depends on the orbital from which the electron is removed. Clearly, in the case of the removal of a valence electron, the relaxation energy is expected to be very small because mainly the electrons above the level from which the electron is removed contribute to screening. On the other hand, for the screening by the cluster valence electrons only the net positive, almost pointlike hole of charge |e| counts, and it should not matter much from which inner-atomic shell a core-electron is removed. Possible differences in the relaxation energy are expected to be due to other core-shell electrons. We call this inner-atomic screening, in contrast to inter-atomic screening by cluster valence electrons. A widely used approach in the DFT community is to use the KS orbital energies as a zeroth-order approximation for the binding energy of electrons. However, the initial state KS orbital energy may be above or below E iσ,rlxd , depending on the xc potential used [11]. Energy differences E are less sensitive and, even more so, binding energy shifts, i.e. the difference of energy differences.

Metal-sphere model
The energy required to remove one electron from a metal sphere of radius R that, initially, has Z e negative elementary excess charges −|e| reads Here, W is the work function. In the bulk limit R → ∞ one has E metsph = W , independent of Z e . Equation (14) can be simply derived by calculating the initial and final energies. As the interior of the metal sphere is, by definition, field-free, one has (with E(r ) being the electric field) so that The work W needed to remove an electron from the bulk is added 'by hand', leading to (14). Improved metal-sphere models yield also W and help to resolve the 'image charge paradox' [21] 2 . Figure 2 summarizes schematically the prediction of the metal-sphere model. Starting from the bulk value W for the binding energy for 1/R → 0, equation (14) predicts a linear behavior in 1/R with a negative slope for initially negatively charged clusters, a positive slope for initially neutral clusters and increasingly steep positive slopes for more and more positively charged clusters. While the cluster radius is-apart from a small spill-out-well-defined for big clusters, this is not the case for small clusters. In the following, we will use (4) down to the single, isolated atom (where R = r s ) and hence plot the respective binding energy versus 1/r s . Even with the uncertainty in the definition of the single atom or ion radius kept in mind, the binding energy for the single atom or ion in general does not lie on the metal-sphere line predicted by (14). Hence, there is a transition region (hatched area in figure 2) in which the binding energy is expected to deviate from the metal-sphere line toward the single atom or ion result (filled colored circles in figure 2). Table 1 summarizes the PZ-SIC xLSD results for the smallest and largest systems we consider in this work, namely the single, isolated Na atom (A = 1) and the Na 58 cluster (A = 58). One sees that in the single-atom case the removal of the 3s↑ valence electron and subsequent relaxation gives E 3s↑,rlxd = 0.191, in good agreement with the experimental ionization potential 0.189. 3 Table 1. PZ-SIC xLSD results for the single, isolated Na atom (left) and the Na 58 cluster (right). The subscripts a and c on the right-hand side indicate whether the level belongs to the central atom or jellium cluster, respectively. The Na 58 cluster is a spin-neutral system with closed shells (i.e. A = 58 is a 'magic number') so that binding energies are independent of whether a spin-up or spin-down electron is removed. The same procedure for a 2p↓ yields a binding energy E 2p↓,rlxd = 1.395, i.e. 38 eV, in very good agreement with the experimental result for the 3 P 2 XPS peak [22] 4 . The removal of a 2p↑-electron gives the slightly different and higher binding energy E 2p↑,rlxd = 1.403. Looking at the KS orbital energies one recognizes that − 3s↑,rlxd = 0.186 is slightly below E 3s↑,rlxd = 0.191. Hence, Koopmans' theorem is almost fulfilled for the valence electron, which is not surprising because relaxation effects are expected to be small if a loosely bound outer electron is removed.

Results and discussion
The difference between KS orbital energy and relaxed energy difference is larger for the 2pelectron removal: − 2p↓ = 1.452 > E 2p↓,rlxd = 1.395. Obviously, the relaxation energy does not fully cancel the 'non-Koopmans energy' [11].
For Na 58 we see on the right of table 1 that the binding energy of the valence electron E 1g,rlxd is reduced compared to the atomic case, as expected from the metal-sphere model for an initially neutral cluster. The binding energy E 2p,rlxd = 1.308 also shifts to lower values with increasing cluster size. However, the absolute value of the KS orbital energy increased with respect to the atomic value. Hence, naively estimating the photoelectron peak position in an XPS spectrum to be kin = 2p +hω is not only inaccurate but may even predict peak shifts into the wrong direction. Red crosses indicate experimental data points, red lines (several dashed in (a) and one solid in (b)) the metal-sphere model prediction (14) for Z e = 0. For the valence case (b), there is almost no difference between relaxed binding energy, unrelaxed binding energy and the KS orbital energy. However, for the removal of a 2p core-shell electron the differences are large and increase with increasing cluster size. The difference between unrelaxed and relaxed binding energy is the relaxation energy, as indicated by the curved arrow in (a). The purple '+' symbols in (a) show the minimum absolute value of the 2p KS orbital energy for the relaxed configuration as the 2p ground state fractional occupation number is varied between 1 and 0 (see the text for discussion). Figure 3(b) shows E valence,rlxd , E valence,unrlxd and | valence | versus 1/R. We observe the wellknown zigzag behavior of the binding energy because of shell-closures [13,23,24] before the metal-sphere result is approached for increasing cluster size. For the removal of a valence electron the KS orbital energies | valence | are very close to E valence,rlxd , and both develop toward the metal-sphere behavior as 1/R decreases. The relaxation energy valence = E valence,unrlxd − E valence,rlxd is very small. The typical valence electron result for binding energies in clusters shown in figure 3(b) has been known for decades [13,23,24]. The qualitatively new aspect coming into play when core-electrons are removed is the reorientation of the AZ + Z e − 1 electrons after the emission of the core-shell photoelectron. There are (at least) three options:

Energies versus 1/R
(i) The photoelectron is slow enough to gain the relaxation energy which is set free due to the screening of the core-hole by the other electrons but fast enough to escape before the core-hole itself decays. (ii) The photoelectron escapes so rapidly that it cannot pick up the relaxation energy. (iii) The photoelectron is slow enough to gain energy not only from relaxation but also from Auger-like core-hole decay.
Signatures in XPS spectra that would support option (iii) were not observed in experiments [16,17]. Instead, slow 'secondary' electrons are measured, which are likely generated via Auger decay of the core-hole after the photoelectron is already gone. Indeed, the estimated photoelectron escape times are sub-femtosecond for the cluster sizes and photon energies considered in this work and thus faster than Auger decay [25].
Option (i) is supported by the experiment because the experimental binding energies (6) for sufficiently big clusters follow again the metal-sphere prediction with only the work function shifted to the respective bulk binding energy.
Option (ii) would manifest itself as an increased binding energy of the photoelectron as compared to the situation where all the relaxation energy goes to the photoelectron. Instead, the experimentally measured binding energies in [16,17] drop below the metal-sphere result, which ultimately must be so because the single-ion value (Pb − in [16,17]) lies below the metal-sphere line. Hence, one has to distinguish (at least) two effects here: the ability of the other electrons to screen and the ability of the photoelectron to pick up the relaxation energy due to screening. The first one seems to be relevant in the experiments [16,17] so that indeed the metal-to-nonmetal transition in Pb − A around A = 20 could be probed. The situation for initially neutral Na clusters is quantitatively different. In fact, we find-apart from the shell-oscillations-little deviation from the metal-sphere behavior down to the single-atom limit in figure 3(b), in agreement with experiment [26]. Figure 3(a) shows the PZ-SIC xLSD result for the removal of a 2p core-shell electron from the central, all-electron Na atom. Again, the relaxed energy difference follows approximately the trend predicted by the metal-sphere model, although with substantial fluctuations on the eV level (1 eV 0.037 au). These fluctuations are also correlated with the filling of electronic shells but in a less obvious way as compared to the valence case where maxima in the (absolute value of the) binding energy occur at half-filled and full shells. It would be desirable to move toward bigger clusters in order to see whether our model really approaches the metal-sphere result for the screened 2p core-hole. However, it is numerically very demanding to proceed significantly towards the bulk limit. Note that 1/R = 0.05 already amounts to A = 125. Figure 3(a) clearly shows the importance of the relaxation energy 2p = E 2p,unrlxd − E 2p,rlxd (i.e. the difference between blue triangles and black squares). As expected, the relaxation energy increases with the cluster size. For the biggest cluster treated in this work, Na 58 , we have (58) 2p = 0.3, i.e. 8.2 eV.
The KS 2p-orbital energies | 2p | of the ground state (green diamonds) are between the relaxed and the unrelaxed energies, as is expected from PZ-SIC xLSD [11]. We checked how the KS orbital energy shifts as the fractional occupation number f 2p is decreased from 1 to 0 in the always relaxed system. We observed that | 2p ( f 2p )| decreases below E 2p,rlxd (purple '+' symbols in figure 3(a)) and then, for 0.2 < f 2p 0, proceeds toward the value for the coreshell-ionized system. If the unknown exact xc potential could be employed the orbital energy would jump discontinuously as f 2p changes from 0 < ε 1 to 0. This is the so-called 'derivative discontinuity' [27]. The PZ-SIC xLSD mimics this behavior but smoothes out the discontinuity. We see that the metal-sphere result in figure 3(a) lies between the relaxed binding energy and the minimum KS 2p-orbital energies min | 2p ( f 2p )|. Fractional occupation numbers are related to Slater's transition state theory where the idea is to consider a (relaxed) configuration 'halfway' between the initial and the final state. Indeed, 2p ( f 2p 0.5) +hω would predict a photoelectron peak in XPS spectra close to what is expected from the metal-sphere model, whereas 2p ( f 2p = 1) +hω would be completely off and-apart from fluctuations-erroneously increasing with increasing cluster size.
Concerning the metal-sphere result for the 2p core-level binding energy, it should be noted that the 2p bulk limit is not precisely known. Moreover, the experimental peaks are typically split over 0.4 eV because of fine structure. We estimated the 2p bulk value from extrapolating the binding energy obtained from the large-cluster XPS spectrum in [22]. As a consequence, the core-level metal-sphere result in figure 3(a) may have an uncertainty of ±1 eV ±0.04 atomic units (red dashed area).
The experimental results for the core-level binding energies as a function of the cluster size [16,17] follow the metal-sphere prediction down to a certain minimum cluster size and then deviate, ultimately approaching the single atom or ion result, which in the experiment with Pb − clusters was below the metal-sphere binding energy for the 5d and 4f core-shells (unlike for the 2p shell in Na above). A too rapid escape of the photoelectron would lead to an increased binding energy. It thus seems that the relaxation energy that is set free due to screening goes indeed to the photoelectron and that the deviation from the metal-sphere model is an elementspecific electronic structure effect. In fact, it is expected that from a certain minimal cluster size on the spherical jellium approximation breaks down. However, our extended jellium model has the exact single-atom limit, so it may deviate at a too small cluster size from the metal-sphere line but it ultimately will. The only exception arises when the single atom or ion binding energy itself lies-perhaps by chance-on the metal-sphere line. Indeed, for sodium this is almost the case, as is seen in figure 3.

Time-dependent density functional theory study of the screening dynamics
In order to understand why the relaxation energy that is set free due to screening goes to the photoelectron, we determined the relevant time scales. To this end, we removed at time t = 0 instantaneously a 2p↓ electron from the central, all-electron Na atom in the Na 58 cluster, and propagated the KS orbitals thereafter. Figure 4 shows the difference in the radial spin densities as a function of time t and radial position r . The spin-down density n ↓ shows a very similar behavior and is therefore not shown. At small r < 1, mainly the charge density of the inner electrons (up to the 2p shell from which the photoelectron has been removed) screens on a  very fast time scale ( 2 au 48 as) and continues to oscillate. This is simple to understand: after the sudden creation of the core-hole, the KS orbitals experience a new potential in which they are no longer eigenstates. If the potential was stationary one would expand the KS orbitals in the new eigenstates, and the oscillation periods would be given by the inverse eigenenergy separations. In TDDFT, the situation is more involved because the potential is nonlinear, i.e. it depends on the KS orbital densities. Moreover, we cannot employ linear response theory here because the change in the fractional occupation f 2p↓ number is not small but unity. Nevertheless, the energy width spanned by the inner electronic levels is such that it yields the screening time scale visible in figure 4. As can be inferred from the density oscillations for bigger r , the valence electrons screen on a time scale of 20 au 0.5 fs, i.e. still sub-fs, in agreement with [28] 5 . The decay of the density oscillations seen in figure 4 is probably underestimated because spontaneous emission, electron ion collisions, energy transfer to ionic degrees of freedom and perhaps other possible dissipation channels are absent in our TDDFT approach.
The time-dependent calculations clearly reveal that the screening dynamics, and thus the relaxation, are fast enough to provide the escaping photoelectron with the relaxation energy. For 5 In [28], the dynamics after the sudden immersion of a negative charge in a jellium cluster was investigated. instance, the kinetic energy of a 2p electron in a Na 58 -cluster that is removed by a 40 eV photon is 4.6 eV, and it takes 0.6 fs to escape from the cluster. Hence, it certainly has enough time to pick up the relaxation energy from its own shell (and the shells below), and even the time scale for the valence electron screening is comparable.

Conclusions
A jellium model with a central all-electron atom has been introduced in order to model recent experiments on core-shell ionization of metal clusters as a function of the cluster size. DFT in xLSD with a PZ-SIC has been employed to calculate the self-consistent electron configurations of the model applied to sodium clusters. The relaxed binding energy follows the trend predicted by the metal-sphere model and confirmed experimentally. However, shell-structure-related fluctuations more pronounced than those for valence electron removal are observed. Both the unrelaxed binding energy and the initial, ground state KS orbital energy are completely off, demonstrating a substantial relaxation energy, the importance of screening and the irrelevance of initial-state core-shell KS eigenvalues. Predicting a photoelectron peak by simply adding the photon energy to the initial KS core-level energy is, in general, doomed to fail. Instead, for the xc potential used in this work KS eigenvalues evaluated for fractional occupation numbers (in the spirit of Slater's transition state theory) are closer to the relaxed binding energy and metal-sphere result.
TDDFT has been applied in order to investigate the screening dynamics after the sudden removal of a Na 2p core-electron. Inner-atomic screening due to electrons in the same and other core-shells of Na occurs within a few tens of attoseconds. Inter-atomic screening by the cluster valence electrons occurs rather on the hundreds-of-attoseconds time scale. The inner-atomic time scale is, of course, core-shell dependent and may be even faster for deeper-lying shells and heavier elements, while the screening time scale in metal clusters due to valence electrons is less sensitive and depends mostly on the number of valence electrons per atom. Whether the photoelectron 'has enough time' to pick up the relaxation energy depends on the ratio of the escape time to the relaxation time. The higher the photon energy and the smaller the cluster the higher should be the net binding energy. However, the deviation from the metal-sphere model observed so far experimentally is most likely dominated by electronic structure effects. In order to verify this assertion, future studies will go beyond the jellium model.