Final-State Simulations of Core-Level Binding Energies at Metal-Organic Hybrid Interfaces: Artifacts Caused by Spurious Collective Electrostatic Effects.

Core-level energies are frequently calculated to explain the X-ray photoelectron spectra of metal-organic hybrid interfaces. The current paper describes how such simulations can be flawed when modeling interfaces between physisorbed organic molecules and metals. The problem occurs when applying periodic boundary conditions to correctly describe extended interfaces and simultaneously considering core hole excitations in the framework of a final-state approach to account for screening effects. Since the core hole is generated in every unit cell, an artificial dipole layer is formed. In this work, we study methane on an Al(100) surface as a deliberately chosen model system for hybrid interfaces to evaluate the impact of this computational artifact. We show that changing the supercell size leads to artificial shifts in the calculated core-level energies that can be well beyond 1 eV for small cells. The same applies to atoms at comparably large distances from the substrate, encountered, for example, in extended, upright-standing adsorbate molecules. We also argue that the calculated work function change due to a core-level excitation can serve as an indication for the occurrence of such an artifact and discuss possible remedies for the problem.

Introduction

X-ray photoelectron spectroscopy (XPS), also known as electron spectroscopy for chemical analyses (ESCA), is a widely used technique to analyze the chemical structure of surfaces and interfaces, providing qualitative and quantitative information about the chemical neighborhood of specific atoms.[1,2] In addition to the chemical environment, core-level binding energies are also influenced by the local electrostatic potential at the position of the excited atom.[3−6] This effect is related to observations for ionic crystals that the differences in Madelung energies between the bulk and the surface can easily amount to ∼1 eV.[7] Electrostatic shifts are also of particular relevance for interfaces in cases where large dipoles occur. This is very common for hybrid organic–inorganic interfaces (relevant, e.g., for the areas of organic and molecular electronics), where interfacial potential shifts are typically associated with collective (also termed cooperative) electrostatic effects.[4,5,8−12] Also in this context, electrostatically triggered core-level shifts on the order of 1 eV have been observed for polar self-assembled monolayers (SAMs) adsorbed on metal substrates.[3,4,13−15] For such systems, the electrostatic shifts can be straightforwardly rationalized by the periodic arrangement of polar entities at the interface. The superposition of their fields causes a step in the electrostatic energy that changes not only the sample work function but also the energetic positions of core levels relative to the Fermi level of the substrate, as described in detail in ref (5).

The interpretation of experimentally acquired spectra frequently relies on first-principle simulations.[2,16] There is a broad range of different approaches to simulate XP spectra.[2,17,18] The most simple strategy inspired by Koopmans’ theorem[19] would be to associate core-level binding energies with the orbital energies of a ground-state calculation.[3,4,6,20−26] This approach, often referred to as the initial-state approach, does not provide absolute values of the core-level binding energies but yields relative shifts between different systems.[2,27−31] This is often sufficient for understanding molecular adsorbates, like self-assembled monolayers with atoms at rather large distances from a metal substrate. There, the trends obtained within the initial-state approach typically agree very favorably with experimental data[3,4,14,20,32] and core hole screening effects can be accounted for by electrostatic models.[33,34] The situation becomes more involved when considering surface core-level shifts, adsorbates in the immediate vicinity of the surface, or the effect of alloying: also there, initial-state calculations have provided valuable insights. Bagus et al.,[2] for instance, showed that surface core-level shifts at Al and Cu surfaces (as representatives for sp and transition metals) are primarily due to initial-state effects.[35,36] The neglect of (potentially site-dependent) screening in initial-state calculations can, however, cause problems: for example, Methfessel et al.[27] found for intermetallic MgAu compounds that the inclusion of screening effects “changes the picture drastically” with a sign change in the shift of the Mg 1s core state; Pehlke and Scheffler[37] observed “remarkably different” screening effects for the two atoms in the surface dimers of Si and Ge; Stierle et al.[38] reported that the inclusion of final-state effects changes the sign of surface core-level shifts in NiAl(110); and Birgersson et al.[39] found that “a large variation exists in the relative importance of initial- and final-state effects for CO on Rh(111)” and stressed the “dangers of interpreting core-level binding-energy shifts in a simple initial-state framework”.

The resulting need to properly account for screening effects has triggered the development of a variety of more sophisticated tools,[17,18] including Green’s function-based approaches,[40−42] techniques based on the GW approximation,[40,43,44] and the explicit simulation of excitation processes,[45] for example, in the framework of time-dependent density functional theory (DFT).[46−49] The most common strategies for complex, extended systems are based (i) either on calculating the difference in total energy between the (fully) ionized interface and the interface in its ground state—typically referred to as Δ self-consistent field (ΔSCF)[17,37,50−54]—(ii) or on the Slater–Janak transition-state approximation.[17,53,55−63] The latter relies on calculating the (Kohn–Sham) orbital energies of a partially ionized core level, whose occupation has been reduced by 0.5 electrons. These approaches are typically referred to as final-state approaches, as they include screening effects of the core hole by the polarization of the electron cloud of the substrate. In passing, we note that this screening can also be accounted for semiquantitatively in initial-state calculations by employing a mirror charge correction proportional to 1/[4ε(z – z0)] (with z – z0 denoting the distance of the excited core level from the image plane of the substrate and ε referring to the dielectric constant of the adsorbate layer).[4,14,20,34,64,65] One would, however, expect that an explicit consideration of screening effects in the quantum-mechanical simulations through final-state approaches should be superior. This suggests that the final-state calculations should typically provide improved results not only for surface core-level shifts but also for adsorbate layers. In the present paper, we will, however, show that the final-state calculations in conjunction with periodic boundary conditions (PBCs) can give rise to a potentially serious complication that is due to the periodic repetition of the core hole in every unit cell. In this context, we will focus on the practically relevant interfaces between metal substrates and physisorbed organic molecules, in which the distance of the generated core hole from the metal surface is comparably large. We will also primarily discuss results obtained within the Slater–Janak transition-state approximation.

Mapping the Interface on a Suitable Model System and the Role of Collective Electrostatics

For understanding possible pitfalls of final-state calculations, it is crucial to consider yet another methodological aspect, namely, how the interface in question is mapped onto an atomistic model system. Usually, one possibility is to apply slab-type calculations, which apply periodic boundary conditions (PBCs) with the metal substrate represented by a two-dimensional (2D) periodic slab consisting of a finite number of metal layers onto which molecules are adsorbed. In the third dimension, periodic replicas of the slab are then decoupled quantum-mechanically by introducing a wide enough vacuum gap and electrostatically by a self-consistently determined dipole layer.[66] Alternatively, the system in question can be modeled by a finite size cluster employing open boundary conditions.

In the following discussion, we will focus on PBC simulations, as they straightforwardly account for collective electrostatic effects, which typically dominate the electronic properties of organic–inorganic hybrid interfaces.[5,8,11,67] These effects arise from the omnipresence of dipoles at surfaces and the fact that a periodic arrangement of dipoles causes a step in the electrostatic energy, shifting the energy landscapes above and below the dipole layer relative to each other (for a tutorial review, see ref (5)). In this context, it should be mentioned that the actual interfacial charge distributions are typically more complex than mere dipoles but the latter term is still used as a “shorthand” for the actual situation, as often the true interface properties can be conveniently mapped onto a dipole model.[5,9] For organic–inorganic hybrid interfaces with typically rather large interfacial dipoles, it is well established that these effects significantly change substrate work functions and the alignment between electronic levels at the interface.[5,8−11] Most important in the present context is that the shift in the electrostatic energy also changes core-level binding energies. This has, for example, been measured as well as modeled for self-assembled monolayers formed by aliphatic[13,15,68,69] and aromatic molecules[3,15,20] containing polar entities embedded into their backbones.

As the above-described effects arise from the superposition of the electric fields of the periodically repeated dipoles at interfaces, their description is intrinsically well compatible with 2D periodic boundary conditions. Thus, it is also not surprising that PBCs are commonly employed when modeling core-level excitations at organic–inorganic hybrid interfaces.[4,5,7,20,61,62,70−73]

In this context, it should be noted that cluster calculations a priori miss the above-described effects, as they do not consider the periodically repeated polar entities present at the interface. As the impact of the interfacial dipoles in neighboring unit cells is primarily electrostatic in nature,[5] it should, however, be rather straightforward to include them via suitable electrostatic embedding schemes[74−78] (in analogy to what is done to account for the Madelung energies when studying ionic crystals[6,7]). Thus, properly corrected cluster calculations might indeed be an interesting strategy for modeling the said interfaces when employing final-state approaches, as will become evident from the data presented in this manuscript.

In fact, the key advantage of PBC calculations described above can become their Achilles heel, when it comes to final-state calculations: while periodic boundary conditions properly capture the impact of periodically repeated polar entities present at the interface, they also periodically replicate the core holes and the corresponding polarization effects in final-state calculations. This then gives rise to artificial collective electrostatic effects. A closely linked problem is that the unit cells in PBC calculations are typically charge neutral to prevent a divergence of the electrostatic energy (although there are also certain approaches to embed charged unit cells within periodic systems[78]). Thus the (half) electron excited from the core level in the final-state calculation is usually not entirely removed from the system but placed into the lowest unoccupied state, which for metal substrates resides right at the Fermi level.[61,62,79−81] Adding a (fractional) charge to a sufficiently thick metal slab should a priori not pose a sizable problem, as it does not noticeably modify the electronic structure of the substrate. Adding the excess charge to an unoccupied state is considerably more serious for semiconducting and insulating substrates, where the additional (half) electron is put into the conduction band (or into unoccupied states of the adsorbed molecules, if they are lower in energy). This is by no means consistent with what is actually happening at the interface as a consequence of the core-level excitation. The problem is further amplified by the fact that at semiconductor surfaces the additional electrons can be artificially localized in diffuse Rydberg orbitals close to the core holes, as described in ref (6). As a remedy to this specific problem, Bagus et al.[6] suggested compensating for the charge of the core hole via the virtual crystal approximation[82−85] rather than adding charge to the conduction band.

This, however, does not solve the problem that for semiconducting as well as metallic substrates the combination of the core hole and the compensation charge creates a polar unit cell, which in conjunction with periodic boundary conditions results in an artificial dipole layer. The resulting spurious collective electrostatic effects can have significant consequences, especially when the separation between the core hole and the compensation charge is large, as will be discussed below. In passing, we note that these are not the only spurious electrostatic interactions PBC calculations can suffer from,[86] as similar complications also occur when modeling charged defects.[87−90]

Methane on Aluminum: A Simple, Yet Instructive Model System for Metal-Organic Hybrid Interfaces

To illustrate the consequences of the spurious dipole layer, we designed a prototypical test system for an interface consisting of a metal substrate and a physisorbed organic molecule. It consists of an Al(100) surface with methane molecules adsorbed in every 2 × 2 surface unit cell (see Figure ). Methane represents the simplest organic molecule. Its frontier levels are far from the Fermi level of the substrate. Thus, Pauli pushback is the only origin of the actual interface dipole[5,91−94] and even when the (half) core hole has been formed, the lowest unoccupied molecular orbital (LUMO) of the ionized adsorbate molecule will (in most cases) remain well above the metal Fermi level. This prevents a transfer of electrons from the metal into the molecular LUMO as a consequence of the core-level excitation. Thus, one avoids Fermi-level pinning,[5,91] which would further complicate the situation (see below). Moreover, the small size of methane confines intramolecular screening effects to a very small volume. Consequently, the simulation artifacts arising from the spurious dipole layer can be shown without significant interference from additional effects, like massive interfacial charge transfer, Fermi-level pinning, or complex dielectric screening. Aluminum has been chosen as the substrate to minimize computational costs for our all-electron calculations. This is important, as to clearly demonstrate the impact of artificial collective electrostatic effects, the largest considered supercell contains 36 adsorbate molecules and 432 substrate atoms. Moreover, the exact nature of the metal substrate should only have a very minor impact on the results obtained here.

Figure 1

Largest supercell of the methane/Al(100) interface investigated in the present study. The black box indicates the smallest considered supercell (2 × 2 surface unit cell of Al(100); base area of 33 Å2), while the colored box shows the biggest supercell (12 × 12, with a base area of 1181 Å2 and containing 36 methane molecules). In the bottom image, only a part of the vacuum gap is shown.

As parent unit cell, we chose a 2 × 2 surface unit cell of Al(100) containing a single methane molecule. For bigger supercells, the parent cell was repeated simultaneously in both the x- and y-directions, creating 4 × 4, 6 × 6, 8 × 8, 10 × 10, and 12 × 12 supercells with identical adsorbate densities and aspect ratios. The surface area of the investigated cells (cf. Figure ) ranges from 33 Å2 (for the 2 × 2 unit cell containing a single methane molecule) to 1181 Å2 (for the 12 × 12 supercell containing 36 methane molecules). We emphasize that the trends discussed below are not dependent on the adsorbate density but rather are determined by the density of (half) core holes created on the surface. This is shown explicitly in the Supporting Information (Section 1), where we compare calculations on supercells containing only a single adsorbate molecule (which is also excited) with simulations on supercells with the same methane coverage as in the parent 2 × 2 unit cell (where only one carbon atom per supercell is excited). Both sets of calculations yield essentially the same results, which implies that the relevant quantity for the calculated core-level shift (at least for the present model interface) is the excitation density rather than the coverage of the adsorbate molecules. In the following, we will report data obtained for the full-coverage case with a single methane molecule per supercell excited, such that the inverse supercell size directly corresponds to the excitation density. The only exception is a 3 × 3 surface unit cell at slightly reduced coverage, which we calculated to generate an additional data point with an excitation density between the 2 × 2 and 4 × 4 cells.

Methodology

The quantum-mechanical modeling was done using density functional theory (DFT) employing the all-electron, full-potential FHI-aims code version 180808.[95−99] To perform the slab-type band-structure calculations on the interfaces, the Perdew–Burke–Ernzerhof (PBE)[100,101] functional was used for describing exchange and correlation. Long-range van der Waals interactions were accounted for by employing the Tkatchenko–Scheffler dispersion correction.[102] All calculations were done with tight settings for all atomic species (as supplied by FHI-aims). Details on the corresponding basis functions are described in the Supporting Information (Section 2.1). Reciprocal space was sampled by a Γ-centered 12 × 12 × 1 grid for the 2 × 2 unit cell and by smaller grids for the larger cells (8 × 8 × 1 for the 3 × 3, 5 × 5 × 1 for the 4 × 4, 4 × 4 × 1 for the 6 × 6 and 8 × 8, and 2 × 2 × 1 for the 10 × 10 and 12 × 12 supercells). These grids are very well converged (despite minor variations in the k-point density between different supercells), as shown in the Supporting Information (Section 2.2). The change in the electron density between subsequent iterations was converged to 10–5 e–, and the change of the total energy of the calculated system was converged to 10–6 eV.

Interfaces were modeled employing the repeated-slab approach with the Al substrate represented by three metal layers. These are rather few layers but making this choice was inevitable to consistently calculate also the largest supercells (the 12 × 12 supercell containing 612 atoms). The Supporting Information (Section 2.3) contains layer-convergence tests for the 2 × 2 unit cell (the cell most strongly affected by the artifacts discussed here), which show that for the current system also three-layer slabs yield reliable core-level binding energies. The periodic replicas of the slabs were quantum-mechanically and electrostatically decoupled by a vacuum gap of 30 Å and a self-consistent dipole correction.[103] The geometries of the adsorbate molecules in the 2 × 2 and 3 × 3 unit cells were fully relaxed until the remaining forces on each atom were below 10–3 eV/Å. The geometries of the other supercells considered here were directly derived from the 2 × 2 cells. The obtained adsorption height amounts to 3.68 Å for the central C atom above the topmost Al layer. Whether the top Al layers are relaxed only negligibly impacts the obtained results (see the Supporting Information, Section 2.4); thus, we stick to unrelaxed surface layers.

The final-state calculations were done within the Slater–Janak transition-state approximation.[55−59] To that aim, half an electron is removed from the carbon 1s core orbital of one of the methane molecules in the supercell. The Slater–Janak theorem has a rigorous theoretical foundation for finite size systems,[104] but when employing periodic boundary conditions, a complication arises, as the unit cell needs to stay charge neutral (see the discussion in Section ). Thus, in that case, the excited charge is moved to the lowest unoccupied level, which for the metal substrate considered here corresponds to a state right at the Fermi level. Notably, the region of electron accumulation following the core-level excitation is found right above the metal surface underneath the excited molecule, as will be discussed below. Correspondingly, one is dealing with one excitation per supercell and an excitation density that is inversely proportional to the size of the supercell. As each core hole excitation creates a dipole at the surface, this inverse relation also applies to the dipole density.

For supercells, special care had to be taken that in the ground-state calculations the orbital from which the charge was eventually removed was localized on only one carbon atom. To achieve that, in most cases, the translational symmetry within the supercells had to be broken by moving one methane molecule by −0.01 Å closer to the surface. The occupation of the selected orbital was then reduced and kept fixed in the following self-consistent field (SCF) cycles. All calculations were performed in a spin-unpolarized manner, as commonly done for such core-level excitation simulations.[6] Spin-polarized calculations on selected test systems yield equivalent trends, as shown in the Supporting Information (Section 2.5).

In addition to employing the Slater–Janak transition-state approximation, we also performed ΔSCF calculations on selected systems, as this approach has also been used repeatedly in conjunction with periodic boundary conditions for calculating core-level excitations at interfaces and surfaces.[17,51,53,54,63] In the ΔSCF approach, the core-level binding energy is obtained as the difference in total energy between the system with, in this case, a full electron excited from the core level and the system in its ground-state configuration. Notably, while in ΔSCF calculations on finite size clusters, the excited electron can be removed from the system, when employing periodic boundary conditions it is again typically put into the lowest unoccupied orbital. If that orbital is localized in the substrate, this again results in a spurious dipole layer. In the ΔSCF calculations, due to the excitation of a full electron, the magnitude of the dipole is essentially doubled compared to the Slater–Janak case (unless this is prevented by pinning effects; see below).

For the analysis and visual representation of the data, Python was used in conjunction with NumPy[105] and matplotlib.[106] Ovito[107] was applied for generating the three-dimensional (3D) view of the system, and VESTA[108] and XCrySDen[109] were used for producing isodensity plots. The figures were compiled with GIMP.[110] The graphs in Figures b and 6 have been compiled using Mathematica, version 11.3 from Wolfram Research.

Figure 2

(a) C 1s core-level binding energy of methane adsorbed on Al(100) as a function of the chosen supercell size for calculations employing the final-state approach within the Slater–Janak transition-state approximation. (b) Electrostatic energy landscape for an electron generated by an isolated pair of a negative and a positive point charge (top panel) and by two oppositely charged, square periodic, 2D arrays of point charges (three lower panels). The distances between the charges in the arrays in the three lower panels scale as 3:2:1, and their packing densities scale as 1/9:1/4:1. While the electron electrostatic energy becomes constant for the isolated dipole in the top panel, there is a step in energy for the pairs of periodic charge arrays. These steps are schematically indicated by the blue arrows.

Figure 6

(a) Dependence of the point-charge-derived correction energy calculated employing eq , Ecorr, on the size of the unit cell, with ε set to 2.1. The vertical lines denote supercells considered in the present manuscript. In the simulation, the actual lattice constant of our model system (5.728 Å), the optimized adsorption distance of 3.678 Å, an image plane of 1.59 Å above the topmost Al layer,[135] and half an elementary charge at every point charge position have been used. (b) Dependence of the point charge-derived correction energy calculated employing eq , Ecorr, on the effective dielectric constant describing screening processes at the interface. The vertical line at a dielectric constant of 2.1 indicates the situation quoted in the main manuscript. The simulations have been performed using Mathematica.[136]

Results and Discussion

Impact of the Size of the Surface Unit Cell on the Calculated Core-Level Shifts

Figure a shows the calculated C 1s core-level binding energies as a function of the supercell size for the methane/Al(100) interface. Applying the Slater–Janak transition-state approximation, they are obtained as the orbital energy of the partially ionized C 1s orbital relative to the system’s Fermi energy. The core-level binding energies vary by as much as 1.2 eV, with the most negative binding energy obtained for the 2 × 2 cell and the least negative binding energy calculated for the 12 × 12 cell. This happens in spite of the identical chemical nature of the studied interface for all supercells. The only appreciable difference between the different supercells is the density of the created core holes (i.e., the excitation density) with only one C 1s core hole generated per supercell.

The screening charge in the metal is found right above the topmost metal layer, as shown in the Supporting Information (Section 3). It is localized exclusively below the ionized methane molecule, such that an effective interfacial dipole is created (which is modified by screening effects within the methane molecule). Due to the required charge neutrality (see above), the screening charge can be associated with the half-electron added to the system right at the Fermi level, although even if the half-electron was entirely removed from the system, the metal would be polarized by the core hole.

The main problem is that the dipole formed by the core hole and the screening/compensation charge in the metal exists in every unit cell due to the periodic boundary conditions. This is in sharp contrast to the actual situation in the experiments, where neighboring core holes with their associated polarization charges are well separated. As a consequence of collective electrostatic effects,[3,4,20] the artificial array of core holes and their countercharges in the metal creates a gradient of the electrostatic energy across the interface. This is shown schematically in Figure b, where we compare the electrostatic energy of an isolated pair of a negative and positive point charge (top panel), with the situation for 2D charge arrays of varying density. While for the single pair, the energy approaches a constant value far away from the point charges, there is a step in the electrostatic energy between the region left of the negative and right of the positive charges for all 2D charge arrays. As a consequence, the array of core holes and polarization charges causes a shift of all electronic states (including the core levels) in the adsorbate layer relative to the Fermi level of the substrate. The magnitude of that shift depends on the density of core holes, as shown in the three lower panels of Figure b. As discussed in the Supporting Information (Section 5), it also scales linearly with the amount of transferred charge (i.e., the effect doubles for a full- compared to a half core-hole calculation). This is insofar relevant, as Williams et al.[57] suggested that removing 2/3 instead of 1/2 of an electron (as in the original formulation by Slater) would yield numerically more accurate values for the core-level binding energies. According to the data in Section 5 of the Supporting Information, this would a priori increase the artificial electrostatic shift of the orbital energy by a factor of 4/3. It would, however, not increase the impact of artificial collective electrostatics on the core-level binding energy, as in the model by Williams et al., the orbital energy for the partially ionized system enters the expression for the ionization energy weighted by a factor of 3/4. Equivalent considerations also apply to approaches relying on the calculation of the slope of the dependence of the orbital energy on the fractional occupation, which have also been suggested by Williams et al. in ref (57).

Notably, not only final-state calculations based on the Slater–Janak transition-state approximation are adversely affected by the presence of the artificial, excitation-induced dipole layer, but also ΔSCF-type final-state calculations are significantly impacted, as is shown in the Supporting Information (Section 6). In that case, the origin of the problem is the additional energy cost associated with the creation of an artificial dipole layer instead of a single dipole.

The above considerations show that both types of final-state calculations (ΔSCF and Slater–Janak) are adversely affected by artificial collective electrostatic effects, which become stronger for high excitation densities. Naturally, the problem can be mended by considering larger supercells. In fact, as shown in Figure , for the system considered here, the core-level binding energy obtained with the 8 × 8 supercell is within ∼0.01 eV of the result for the 10 × 10 and 12 × 12 cells. This means that for a methane molecule in the immediate vicinity of the metal substrate, the calculation of the 8 × 8 supercell containing 16 molecules can be considered to be energetically converged.

Impact of the Creation of the Core Hole on the Global Energy Landscape

As a next step, we discuss the direct impact of the excitation-induced dipoles on the electrostatic energy in more detail. As a starting point for that discussion, Figure a shows the plane-averaged electrostatic energies for the ground-state configurations of the 2 × 2 and 12 × 12 supercells. They coincide, underlining the identical physical and chemical natures of the two systems. The minima in the electrostatic energy at the positions of the Al planes as well as in the region of the adsorbate molecule are well resolved, and the minor methane-induced shift in the sample work function due to Pauli pushback[5,91,111,112] can be inferred from the energetic difference of the vacuum levels left and right of the slab (indicated by the dashed black line).

Figure 3

(a) Plane-averaged electrostatic energy of the methane/Al(100) interface for the smallest (2 × 2; dotted blue line) and largest (12 × 12; solid pink line) cells. The work function on the methane side of the slab amounts to 4.31 eV, while on the Al side it is 4.40 eV (yielding an adsorption-induced work function change of 0.09 eV; see the dashed horizontal line). (b) Calculated change in the plane-averaged electrostatic energy between a final-state calculation (including a half core-hole excitation), Ees,fs, and the ground state, Ees,gs for different excitation densities caused by different supercell sizes. The latter are denoted directly in the graph (the lines for the 8 × 8 (violet) and 10 × 10 (brown) cells are not labeled due to the limited available space). (c) Work function, Φ, of the systems with a half core-hole excitation per supercell on the methane side of the slab. The work function obtained in the ground-state calculation is indicated by the dashed black line at 4.31 eV. The work functions on the Al side essentially do not change, as shown in the Supporting Information (Section 7). The values obtained with the simple electrostatic model in the main text are indicated as short horizontal lines for the corresponding supercells.

The situation changes fundamentally in the presence of core hole excitations (again half a core hole per supercell). This is shown in Figure b by the modification of the electrostatic energy due to the excitation of a half core-hole per unit cell (including the corresponding compensation charge in the metal). In particular, one sees that for smaller supercells (i.e., higher excitation and dipole densities in the final-state calculations) the change in the averaged electrostatic energy is more pronounced. Notably, the majority of that change occurs between the region of electron accumulation right above the topmost metal layer and the position of the carbon atom of the adsorbed methane molecules. When analyzing the work function of the system in the presence of a core hole excitation (given as the difference between the Fermi energy and the electrostatic energy far above the surface), one observes a strong, supercell-dependent change, as shown in Figure c. Even though this work function is not of immediate relevance for the studied interface, it is still useful for judging the collective electrostatic artifact: the overall trend for the evolution of the work function shift in Figure c is similar to that for the core-level binding energies in Figure a and, as discussed in the Supporting Information (Section 4), it is essentially inversely proportional to the size of the unit cell (depolarization effects notwithstanding).

In fact, it is even rather straightforward to obtain a rough estimate of the work function change due to the core-level excitation solely based on electrostatic and geometric arguments, i.e., without performing any quantum-mechanical simulations (see the Supporting Information, Section 4.2). Such a model yields the values indicated by the short horizontal lines in Figure c, which agree rather well with the actual work function changes and, thus, provide a first handle for estimating the adverse impact of the artificial electrostatic effects discussed here.

In this context, it is, however, worthwhile mentioning that the absolute magnitude of the work function shift is significantly larger than the shift of the core-level binding energies (as can be inferred from a comparison of Figures a and 3c). To understand that, one has to analyze the difference between the electrostatic energy at the location of the core level and in the far-field, high above the interface.

Local Impact of the Dipoles

As a first step for analyzing locality effects, it is useful to provide a spatially resolved illustration of the differences in the change in electrostatic energy that occurs due to the presence of an isolated core hole and due to a dense layer of core holes. To that aim, Figure compares the 2D cross sections of that energy change for the 2 × 2 and 12 × 12 cells. While for the 12 × 12 supercell, one essentially obtains the situation of an isolated dipole, which affects the energy landscape only locally, for the 2 × 2 cell, the above-discussed global shift of the electrostatic energy also (infinitely) far above the metal surface occurs.

Figure 4

Calculated difference in electrostatic energy between the final-state and the ground-state calculations for the 2 × 2 (top) and the 12 × 12 (bottom) cells. The values are plotted for a plane parallel to one of the unit cell axes, containing the nuclei of the C atoms. The overlaid atomic structure of the interface is shaded in the right half of the plots to better resolve changes in electrostatic energy close to the nuclei. Furthermore, only part of the 12 × 12 supercell is shown and the 2 × 2 unit cell is repeated four times in the vertical direction.

Figure also reveals that the change in electrostatic energy varies significantly in the lateral direction, especially at low coverage. A detailed analysis of the impact of these lateral fluctuations shows that they are responsible for the much smaller artificial shifts in core-level binding energies compared to the shifts in work functions (see above). This is a direct consequence of core-level binding energies being sensitive to the local electrostatic energy at the position of the orbital from which the electron is excited.[3,8,20] In contrast, the work function measures the electrostatic situation in the far-field, such that the lateral fluctuations of the potential energy are no longer relevant. This aspect is of practical relevance as it, for example, results in core-level shifts in polar SAMs becoming sensitive to sample inhomogenieties.[20] In the present context, it means that averaging the electrostatic energy over the entire unit cell (like in Figure b) very well characterizes the situation for the work function but fails to properly capture the “locality” of core-level shifts. To show that also quantitatively, one has to analyze the electrostatic energy averaged over a much smaller area than in Figure b. This is done in the Supporting Information (Section 8), with the results fully supporting locality as the main reason for the smaller artifacts in the core-level energy calculations.

Impact of the Adsorption Distance on Core-Level Energies

The next aspect to be addressed is to what extent the discussed artifact depends on the geometrical structure of the interface. In particular, we will address how spurious collective electrostatic effects depend on the distance between the metal surface and the atom probed by XPS. This is, for example, relevant when bulky side groups lift the backbone of an adsorbed molecule from the metal substrate. It is also relevant for upright-standing adsorbates,[113−119] e.g., in self-assembled, covalently bonded monolayers.[120−126] To qualitatively assess the situation, we systematically varied the distance between the methane molecule and the surface. In this way, effects arising from (system-specific) screening by the upright molecular backbones between the probed atom and the substrate are not accounted for. Still, the model serves to illustrate the impact of an increased charge-transfer distance.

Figure shows the dependence of the core-level binding energy on the adsorption distance of the methane molecules (specified relative to the equilibrium distance) for differently sized supercells. This plot reveals several interesting aspects:

Figure 5

(a) C 1s core-level energies for methane on Al(100) relative to the Fermi level for different unit cell sizes and different adsorption heights of the methane molecule (blue: 2 × 2 unit cell, orange: 3 × 3 supercell, green: 4 × 4 supercell, red: 6 × 6 supercell, violet: 8 × 8 supercell, brown: 10 × 10 supercell). No data points for 12 × 12 supercells could be included in this plot, as for this very large cell we failed to converge the SCF cycle in the calculations for increased adsorption distances. (b) Difference in core-level energies between the 10 × 10 and the 4 × 4 supercells.

The data points for the 10 × 10 and 8 × 8 supercells essentially coincide at small distances. This implies that there the 8 × 8 and 10 × 10 data display the true effect of core hole and screening. For larger distances, the deviations between the 10 × 10 and 8 × 8 unit cells increase, and at 3 Å above the equilibrium distance (i.e., at an adsorption distance of 6.69 Å) they exceed 0.1 eV; notably, at that adsorption distance also the excitation-induced work function change reaches 0.5 eV, even in the 10 × 10 cell. This is attributed to the growth of the dipole per unit cell upon increasing the charge-transfer distances. In view of the observation that at large adsorption distances even the 10 × 10 data are impacted by artificial collective electrostatic shifts, we refrain from fitting the data with simple electrostatic image charge corrections (where the shift of the core-level binding energies would be inversely proportional to the distance of the excited atom from the mirror charge plane).[34]

The calculated distance dependence of the core-level binding energies clearly deteriorates for the 6 × 6 and 4 × 4 unit cells. In fact, the deviations between the 10 × 10 simulations and the 4 × 4 simulations skyrocket for larger distances, as is shown in Figure b. These results imply that to mitigate the adverse impact of artificial collective electrostatic effects, at larger adsorption distances, one would have to simulate even larger supercells than those studied here. Considering that already the calculation of the 10 × 10 cell is reaching the limits of present computational capacities, studying further enlarged cells is far from trivial. In fact, we failed to achieve convergence in the SCF procedure for 12 × 12 cells at larger distances.

For the smallest considered cells, the situation changes fundamentally. Especially for the 2 × 2 unit cell, the core-level binding energy becomes essentially independent of the adsorption distance. This behavior is reminiscent of the situation encountered for Fermi-level pinning,[5,30,91,112,127−132] as for the smallest unit cells, the shift in electrostatic energy due to the artificial dipole layer becomes so strong that the lowest unoccupied states of the adsorbed molecules would get pushed below the Fermi level of the substrate. This is counteracted by electrons being transferred from the metal into the formerly unoccupied molecular states, significantly modifying the net charge rearrangements (see the Supporting Information, Section 9). The ensuing counterdipole prevents any further increase of the potential step across the interface. Consequently, the energetic positions of the electronic states become independent of the adsorption distance. Such a behavior is inconsistent with the experimental situation for two reasons: first, at least in the present model system, the pinning situation is solely a consequence of the artificial collective electrostatic effects resulting from unrealistically high excitation densities in the simulations. Second, even if the core hole locally shifted unoccupied states below the Fermi level (for example, because of a much smaller energy gap of the adsorbate), for many interfaces, the time scales of the photoelectron experiments[133,134] would be such that the (partial) filling of these states with electrons would be too slow to affect the measured kinetic energies of the escaping electrons (and, thus, the core-level binding energies).

An aspect that is somewhat surprising about the pinned situation for the 2 × 2 unit cell is that pinning occurs at core-level binding energies that are significantly less negative than for some of the considered 3 × 3 and 4 × 4 cells. On the one hand, we attribute this to the different degrees of localization of the core levels and the much more delocalized unoccupied states at which pinning occurs. Consequently, the different orbitals “probe” shifts in the electrostatic energy in different spatial regions and are, thus, differently affected. On the other hand, we observe a more broadened low-lying unoccupied density of states in the 2 × 2 cell compared to the larger supercells. This can also result in Fermi-level pinning already for smaller energetic shifts. The relevant factor here is that pinning occurs at the unoccupied states primarily localized at the excited molecules, as these lie lowest in energy due to the dipole-induced shifts. For the 2 × 2 cell, all adsorbate molecules are excited and, thus, all unoccupied states are shifted by the same amount, favoring a strong intermolecular coupling. Conversely, for 4 × 4 and larger supercells, all excited molecules are surrounded by molecules in their ground state, effectively preventing such a coupling. A more detailed discussion of the impact of the excitation density and the adsorption distance on the shape of the unoccupied density of states can be found in the Supporting Information (Section 9). There, one also finds a discussion of the impact of the choice of the basis set on the unoccupied states of methane, which can quantitatively (albeit not qualitatively) change the situation.

Possible Strategies for Avoiding Spurious Electrostatic Effects in Final-State Calculations Employing Periodic Boundary Conditions

Considering that the discussed artifacts are electrostatic in nature, the question arises whether one could also devise an electrostatic correction scheme. The simplest approach would be to employ a plate-capacitor model like that used in Section . Such a correction only describes the situation in the far-field (i.e., at a sufficient distance from the interface) and does not reproduce the shift of the local electrostatic potential at the position of the excited atom. The latter is, however, what counts for the core-level shifts, as discussed in Section .

A solution to the locality problem would be to explicitly consider the spurious shifts in electrostatic energy due to the periodic replicas of the (half) core holes and their image charges.[87,103] In view of the employed dipole correction (see Section ), periodicity is considered here only in the directions parallel to the surface. The resulting correction for the electrostatic energy of an electron at the position of the central core hole placed at the origin of the coordinate system (at i = j = 0), Ecorr, is then given byHere, a square lattice (like in the studied model system) with lattice constant a is assumed. For the half core-hole calculations, the charge of the core hole and the mirror charge are both set to Q = 0.5e (with e being the elementary charge). D corresponds to the distance between a core hole and its mirror charge and is given by D = 2(z – z0), where z denotes the position of the core hole above the top metal layer and z0 refers to the position of the image plane set to 1.59 Å for Al(100).[135] ε is an effective dielectric constant that describes the screening within the adsorbate layer. When setting ε to 2.1, one obtains correction energies of 1.20 and 0.17 eV for the 2 × 2 and 4 × 4 supercells, respectively. This, indeed, provides an excellent correction of the artificial shifts in Figure (see also Figure a). A complication in this context is, however, that the correct value of ε is a priori not known, while especially for small unit cells its value rather significantly impacts the correction (see Figure b).

This raises the question of whether one could avoid the spurious dipoles in neighboring unit cells by achieving charge neutrality in every unit cell via some static compensation charge rather than by placing the excited (half) electron into an unoccupied level. A homogeneous background charge often employed in bulk calculations[137] will typically still create an artificial dipole layer, whose magnitude depends on the size of the unit cell and, thus, on the size of the vacuum gap separating periodic replicas of the slab.[83] This should become particularly problematic when comparing core-level excitations from atoms at different positions within the unit cell. Another possibility would be to account for the compensation charge via the virtual crystal approximation.[82−85] As shown by Bagus et al.,[2] this eliminates certain artifacts occurring when the compensation charge is put into the conduction band of a semiconductor substrate (see Section 1). As it does not eliminate the spurious interfacial dipole layer, it, however, also does not remove the spurious potential drop between the substrate and the adsorbate.

A strategy for avoiding that drop would be to put compensation charges into a charged sheet above the adsorbate layer (in analogy to the charge reservoir electrostatic sheet technique (CREST) used to model the band bending in the extended depletion regions of semiconductor surfaces; albeit there, the charged sheet is placed below the slab).[138,139] This would, however, also result in an incorrect description of the screening potential between the actually generated core hole and the metal substrate, as the charged sheet tends to concentrate the electric field in the region between the core hole and the sheet. In the region of the metal, this results in a too small field due to the core hole, with the deviation being particularly large for smaller unit cells (like the 2 × 2 cell, as shown in the Supporting Information, Section 11).

Likewise, the problem could not be solved by putting the compensation charge into some higher-lying unoccupied levels within the molecules. This would largely avoid the artificial dipoles, but it would also render all molecules essentially charge neutral (including the one that is actually excited); i.e., in such a calculation, screening by the metal would be eliminated altogether, which would again not represent the actual situation encountered in an experiment.

Thus, we envision four strategies (with varying strengths and weaknesses) that have the potential to either avoid or correct for the spurious collective electrostatic effects discussed in this paper:

Employing the electrostatic correction scheme described at the beginning of this section.

Increasing the supercell size in final-state simulations applying periodic boundary conditions until the artificial core-level shifts are below a certain threshold (see above). In this process, the coverage has to be kept at its initial value while exciting only one atom per supercell. Converging the supercell size will typically be feasible for flat-lying adsorbates, where charge-transfer distances (and, thus, spurious dipoles) are rather small (albeit at significantly increased computational costs). For studying upright-standing adsorbate layers of spatially extended molecules, which are often relevant for practical applications, the required supercells will, however, often be computationally not accessible.

Especially for atoms far from the surface, initial-state calculations with an a posteriori mirror charge correction proportional to 1/(4ε(z – z0)) are a viable strategy to account for screening effects. We have, in fact, routinely and successfully employed this approach for modeling core-level shifts in self-assembled monolayers.[4,14,20,32,68] This approach, however, poses the disadvantage that it describes screening on purely electrostatic grounds, disregarding, e.g., quantum-mechanical effects, which might become relevant especially for atoms that are part of or close to the surface (see Section 1). Moreover, the determination of the relevant parameters (like the dielectric constant of the monolayer) is not necessarily straightforward.

Another strategy to avoid spurious electrostatics would be to refrain from employing periodic boundary conditions and to study finite size clusters as a model for the interface. True collective electrostatic effects at the interface could then be accounted for by properly designed electrostatic embedding schemes,[76−78] like the “periodic electrostatic embedded cluster method” (PEECM).[74,75] As far as the properties of the extended interface vs the cluster are concerned, such corrections would, however, only account for purely electrostatic effects.

Conclusions

The above data show that spurious electrostatic effects can come into play when simulating XP spectra within the final-state framework in combination with periodic boundary conditions. The screening/compensation charge in the metal creates a dipole, which due to the periodic boundary conditions is spuriously repeated in every unit cell. This creates an artificial, dense layer of dipoles, which shifts all electronic states in the adsorbed molecule relative to the metal Fermi level. Consequently, also their core-level binding energies are increased. Therefore, one observes a pronounced shift in the calculated core-level binding energies as a function of the employed supercell size (and the resulting excitation density). For the methane/Al(100) system studied here, this shift amounts to 1.2 eV when comparing calculations for 2 × 2 cells (with an area of 33 Å2) and 12 × 12 cells (with an area of 1081 Å2), when in each cell a single molecule is ionized. A similar effect is obtained when comparing the total energies of the systems in their ground and excited states (i.e., when applying the ΔSCF procedure; see the Supporting Information, Section 6).

Another quantity that is significantly changed by the interfacial dipole layer due to the core-level excitations is the work function obtained in the final-state calculations. Even though this quantity is not of practical relevance for the actual experimental situation, it is still useful for illustrating the artifacts discussed in this paper. Moreover, a significant work function change in a final-state calculation is a strong indication for artificial core-level shifts. This has the advantage that one can predict, whether problems occur, already on the basis of a single calculation without the need to converge supercell sizes.

When increasing the distance of the excited core hole from the metal substrate, the artificial collective electrostatic effects increase even further, such that performing trustworthy final-state calculations in conjunction with periodic boundary conditions becomes virtually impossible.

For cases in which converging the lateral size of the unit cell in final-state calculations becomes impossible, we suggest a simple electrostatic correction, whose application, however, requires a reasonable guess for the effective dielectric constant of the interface. Alternative strategies for studying such situations would be to perform initial-state calculations with an a posteriori mirror charge screening or to join cluster calculations with electrostatic embedding schemes, where each of the suggested approaches has its own strengths and limitations.