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.


Dependence of the core-level binding energies on molecular coverage
In the main manuscript we show that the C 1s core-level binding energies of methane adsorbed on a Al(100) surface calculated employing the final-statefinal-state approach strongly depends on the density of core holes present in the simulation. These results (shown in Figure 2a of the main manuscript) have been obtained when modeling the binding energies as a function of supercell size with one excitation per supercell. In these simulations, the coverage of the methane molecules was kept constant. In Figure S1 we show that essentially the same results are obtained when also varying the coverage of the methane molecules (i.e., when changing the methane density for different unit cells -see Figure S2). The bigger unit cells with different coverages were created by multiplying the 2×2 unit cell and removing every forth methane molecule in the system for the ones with a coverage of 25% and by removing all but one molecule in the case of the 10×10 system, creating a coverage of only 4.0%. For the 3×3 unit cell an aluminum layer with 9 surface atoms per unit cell was created. The adsorption site of the methane molecule was obtained in the same way as in the full coverage 2×2 case, namely by placing the methane molecule on the surface and fully relaxing its geometry (keeping only the positions of the Al atoms fixed). The resulting coverages are summarized in Table S1. This shows that what primarily counts for the described effect is the density of excited core holes.

S4
Figure S1. C 1s core-level binding energy of methane adsorbed on Al(100) as a function of the chosen unit cell size with different methane coverages (see Table S1). The calculations were done employing the final-state approach within the Slater-Janak transition-state approximation. Figure S2. Left

Information on the employed basis functions
The basis functions employed in the FHI-aims simulations have the format in spherical coordinates (r, Θ, Φ) relative to a given atomic center. FHI-aims provides for every atomic species a preconstructed species_defaults file. The used tight basis sets were not further adjusted, because they afforded the required accuracy and efficiency. Note: If a higher tier for the basis set, i.e. when using tight settings, is used, all lower basis functions must be used as well.

Convergence of the k-point grid
The Γ-centered k-point mesh was evenly split along the reciprocal lattice vectors of the unit cell with the same number of k-points in x-and y-direction and 1 k-point in zdirection. This was done because of the quadratic unit cell and the repeated slab approach.
When applying the repeated slab approach, only 1 k-point in z-direction is used, as the interface is not periodic in that direction. Due to the fact that the k-point mesh samples reciprocal space, smaller unit cells require the use of more k-points. Therefore, when the unit cell size was doubled in a given direction, the number of k-points in that direction was halved. To find the number of k-points in x-and y-direction per unit cell length, which are needed to get converged results, several calculations were done with an increasing a As described in the FHI-aims manual, version January 23, 2017.

S7
number of k-points for the two smallest systems until the obtained values of interest did not change beyond a certain threshold. This procedure yielded that 400 k-points per topmost aluminum atom -i.e. 4 Al atoms at the surface in the case of the 2×2 unit cell, 16 surface atoms in the case of the 4×4 supercell, and 9 surface atoms in the case of the 3×3 unit cell -were enough to get orbital energies converged to at least ±0.002 eV. The exact number of k-points used for each supercell and the k-point density per surface atom is shown in Table S3.

Impact of the number of layers contained in the slab
Owing to the very extended supercells required for the present manuscript, all data presented in the main manuscript rely on metal slabs consisting of only three Al layers.
To ensure that this is sufficient for calculating the core-level binding energies of the present system (where there is no significant substrate/adsorbate charge transfer) test calculations for the smallest 2×2 unit cell containing 6 Al layers were also performed. As can be seen in Table S4, the C 1s orbital energies for the two slab thicknesses differed by only ~ 0.01-0.02 eV, even in cases, where the methane atom was moved farther away from the surface slab resulting in a bigger dipole. This suggests that the consideration of threelayer slabs is sufficient for the present purpose. S8

Impact of the geometric relaxation of the topmost metal layer
The aluminum slab was constructed with ASE 2 , utilizing the function to create a fcc (100) surface slab with a tabulated lattice constant of 4.05 Å. To mimic the metal slab, the lowest layer of aluminum atoms was always fixed. To test whether the other layers need to be relaxed, two additional systems were calculated. One with the topmost layer and one with the two topmost layers allowed to relax during a full geometric optimization. Again, this hardly affected the core-level binding energies (see Table S5). Thus, for all calculations the slab with the aluminum atoms in the ideal bulk lattice positions was used.

Spin polarized calculations
To test whether one needs to perform spin polarized calculations, selected supercells were calculated with spin unrestricted settings. Therefore, in the control.in file to start the FHI-aims 1 calculations, the spin parameter was set to collinear and the force_occupation parameter was set by the following lines in the case of a calculation of the 4×4 system: Additionally, an initial guess for the spin of the excited atom needs to be set with the keyword initial_moment in the geometry.in file in the line following the specification of the excited carbon atom. The results obtained with unrestricted spin settings yielded the same qualitative trend (see Figure S3) as the spin restricted calculations. In fact, the mean value of the spin up and spin down orbitals is close to the energy of the C 1s orbital in the restricted calculations.
What is more relevant in the present context is, however, that the energy difference between the 2×2 and 4×4 unit cells is essentially the same for the spin up, spin down, and spin unrestricted channels (1.017 eV; 1.019 eV; 1.019 eV).

Charge rearrangements due to the excitation
When utilizing the Slater-Janak transition-state theory, half an electron is removed from the excited core hole. This charge is then placed into the lowest unoccupied orbital in the calculated system. In the case of a system with a metal substrate, as considered in this manuscript, this 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 shown in Figure S4. Please note that the blue feature in Figure S4 in the region of methane is a consequence of electron depletion due to the polarization of the molecule. Consequently, there is one excitation per supercell and the excitation density is inversely proportional to the size of the chosen supercell. S11 Figure S4. Isodensity plot -front view (a) and top view (b) of the 12×12 supercell showing the electron depletion (blue) and accumulation (red) due to the removal of half an electron from the carbon 1s orbital of the excited molecule. This charge is moved to the lowest unoccupied level in the system by FHI-aims. As can be seen this is localized directly underneath the excited molecule.

Plane-integrated charge rearrangements
Notably, the shapes of the calculated plane-integrated charge rearrangements are independent of the supercell size (cf., data for the 4×4 and 12×12 cells), with minor deviations for the smallest 2×2 cell (solid blue line in Figure S5). The latter we attribute to a somewhat different electronic structure in the 2×2 case, in which all adsorbate molecules are excited. In contrast, in the calculations of larger cells the excited methane molecule is surrounded by inequivalent molecules (i.e., by molecules in their electronic ground state), which minimizes hybridization effects between the molecules. This aspect will be discussed in more detail later, when discussing Fermi-level pinning effects.

Plane integrated rearrangements as a function of the adsorption distance
In the previous section we show the charge rearrangements (Δρfs-gs) in final-state calculations relative to the ground state charge density for different excitation densities (cf. Figure 5). Here we show the charge rearrangements due to different adsorption S13 heights. As can be seen in Figure S6 the charge density difference in and right above the metal slab stays largely the same, no matter how far away the excited methane molecules are moved away from the substrate. Furthermore, the qualitative (and quantitative) charge rearrangements at the molecule prevail, when it is moved away from the metal slab.  Figure S7 shows the integral of Δρfs-gs (the quantity plotted in Figure S6) over the z coordinate of the unit cell as a function of the position up to which the integration has S14 been performed. This integral can either be interpreted as a quantity proportional to the electric field resulting from Δρfs-gs, or, more useful in the present context, as ΔQfs-gs, the cumulative charge rearrangement, which quantifies, how many electrons have been transferred from below to above a plane at position z due to the excitation. Figure S7 shows that the maximum cumulative charge rearrangement somewhat increases with increasing adsorption height, which we attribute to the fact that at larger adsorption distances the overlap between the screening charges in the metal and the polarization charges in the adsorbate diminishes. This observation explains, the slightly superlinear increase of ΔEC 1s with adsorption distance observed in Figure 6b.

Shift of electronic levels
Due to collective electrostatic effects arising from the artificial dipole array in the calculations, all electronic states in the adsorbate layer are shifted relative to the Fermi level of the substrate.

Work function shifts
As a result of collective electrostatic effects, also the work function is shifted depending on the supercell size, i.e. the excitation density. The impact of this shift on the work function is inversely proportional to the base area of the supercell. In Figure S8 the work function change due to the artificial dipole layer in the final-state calculation is plotted over the inverse supercell size. As one can see all datapoints are essentially linear with the exception of the 2×2 unit cell, for which the evolution becomes sublinear. This is due to depolarization effects 3-8 . S16 Figure S8. In the case of Slater-Janak transition-state calculations, Q is set to 0.5 times the elementary charge. For the system discussed in the present manuscript, z amounts to 2 Å as the approximate distance over which the potential drops in Fig. 4b, and r is set to 1.5 (to account for screening effects within the methane layer).

Carbon 1s core-level binding energies
Furthermore, also the core-level binding energies are shifted. The magnitude of that shift again depends on the dipole density, i.e. the excitation density, which is inversely proportional to the size of the supercell. In Figure S9 the core-level energies are plotted over the inverse supercell size. As one can see, in this case an essentially linear relation is obtained. The "linearity" condition is, however, less well fulfilled than for the workfunction change. We attribute this to the "locality" character of the core-level shifts, i.e., S17 the fact that for core-level binding energies also lateral variations in the electrostatic energy are highly relevant.  Figure S11 compares the carbon 1s orbital energies calculated with two different finalstate methodologies, namely the Slater-Janak transition-state theory and the ΔSCF method. In the former, half an electron is removed from the core hole and placed in the lowest unoccupied orbital of the system (in this case at the Fermi level, i.e. somewhere at the metal's surface) and in the latter a full electron is removed from the core orbital and in our case placed into the lowest unoccupied orbital. In the ΔSCF method, the core-level binding energies are then associated with the energy differences of the systems in their S19 excited and ground state, respectively. For all but the smallest unit cell, the core-level binding energies obtained with the two different methodologies are essentially rigidly shifted relative to each other. This shows that the artificial collective electrostatic shifts are similar for both methods. This is insofar intrieguing, as the artefacts play out differently in the two approaches: When applying the Slater-Janak transition-state theory, the collective electrostatic shifts directly affect the energetic positions of the C 1s orbitals relative to the Fermi energy. Conversely, in the ΔSCF method they manifest themselves as the "charging energy" of an interfacial capacitor (cf. main manuscript).

Final-state calculations relying on the ΔSCF approach
The only exception from the overall trend is the 2×2 unit cell, for which the C 1s core-level binding energy calculated within the ΔSCF approach is shifted significantly less than when applying the Slater-Janak transition-state theory. The reason for that is Fermi-level pinning similar to the effects observed in Figure 5a of the main manuscript. The reason why Fermi-level pinning for the ΔSCF approach occurs alrady at the equilibrium distance, while in Figure 5a it occurs especially for larger distances is the increased amount of charge transfer in the ΔSCF (a full electron instead of half an electron), which results in an increased dipole density.  is affected by the added charge. This is done for different supercell sizes (and, thus, excitation densities). The results shown in Table S6 confirm that in none of the considered cases the bottom-side work function is changed significantly, further confirming the notion that the 3-layer Al slab is sufficiently thick for the studies reported in the present manuscript.

Impact of locality effects of the electrostatic energy on corelevel binding energies
Here we show quantitatively that for properly capturing the impact of collective electrostatics on core-level shifts, one has to average the electrostatic energy over much smaller areas than the cross-section of the unit cell. Correspondingly, Figure S12 shows the electrostatic energy averaged over an area of only 0.08 Å 2 parallel to the surface (as an estimate for the "extent" of a C 1s orbital; see section 8.1) for the 2×2 and 12×12 cells.
Both graphs reveal the pronounced dip in electrostatic energy around the core hole. Far enough above the methane molecule, the averaged electrostatic energy gradually approach the far-field values, like in Figure 3b, albeit at a much slower rate. The latter is fully consistent with the data in Figure 4. A disadvantage of the way the data are plotted in Figure S12a is, however, that the massive drops in the energy around the carbon atoms obscures the differences between different excitation densities.
Thus, in Figure S12b we plot the difference between the two curves contained in Figure   S12a. Then, the 12×12 supercell can be viewed as the situation in the absence of artificial collective electrostatic effects (due to the highly diluted dipole density), while these effects are maximized for the 2×2 cell. Consequently, this plot illustrates their impact on a lateral length-scale consistent with the extent of the C 1s electron. In this way it accounts for the local nature at which core-level excitations probe the electrostatic landscape. From these data one learns that the artificial shift in electrostatic energy at the position of the C atom is, indeed, way smaller than the shift of the energy far above the interface. This reconciles the data from Figures 2a and 3c. In particular, it confirms the notion that the different magnitudes of the artificial shifts of the core-level binding energies (in Figure 2) and the work-function changes (in Figure 3c) solely arise from different "locality" with which these quantities probe the electrostatic energy across the interface.

Dependence of evolution of the local electrostatic energy on the size over which the lateral averaging occurs
To approximate electronic wavefunctions, Slater proposed the use of hydrogen-like orbitals with an effective nuclear charge. He also compiled rules for aquiring the effective nuclear charge (Zeff) that should be used to take the screening of the core electrons into account. For a carbon 1s orbital this effective nuclear charge is Zeff = 5.7, 9 which results in a radius of ~0.16 Å within which 2/3 of the C 1s charge density are contained. From this an effective cross-sectional area of the C 1s orbital of 0.08 Å² can be estimated analogous to the area used in the main manuscript (although there, for technical reasons, the averaging occurs over square voxels). In fact, when calculating the plane-averaged charge density associated with the C 1s orbital in the actual interface (see Figure S13) a consistent value is obtained. Figure S13. Plane averaged charge density of the C 1s orbital. The vertical grey lines indicate the extent of the voxels of the grid on which the charge density was written out.

S24
The exact size of the area over which the electrostatic energy is averaged is, however, only of minor relevance: As shown in Figure S14 the depth of the potential well associated with the core hole, as expected, depends strongly on the area over which the averaging occurs.
The actually relevant quantity, namely the difference in averaged potentials for the 2×2 and 12×12 supercells shown in Figure S15 is, however, not affected significantly.

Fermi-level pinning for large adsorption heights and small unit cells
As can be seen in Figure 5a of the main manuscript, for the 2×2 unit cell, the adsorption distance has hardly any impact on the C 1s energies. Conversely, for the 3×3 unit cell upon increasing the adsorption distance the orbital energies are shifted up to a certain value (significantly more negative than for the 2×2) and then pin. In other words, the pinning S26 occurs not only at different distances, but also at very different energies. I.e., a priori, it is not clear, why the core-level in the 3×3 unit cell and in larger supercells can get shifted to energies that are significantly more negative than the energy at which pinning already occurs in the 2×2 unit cell. We here propose differences in orbital hybridization as the explanation for this effect. Figure S16 shows the occupied and unoccupied density of states (DOS) for the 2×2 unit cell at the equilibrium adsorption distance for the ground state and for the final-state calculations. Here we also plot the unoccupied density of states at energies higher than the vacuum level above the adsorbate layer (indicated by colored vertical lines in the following plots). A priori the corresponding states cannot be properly described by an orbital-centered basis set. Nevertheless, we include it here, as it still allows the assessment of similarities between different systems (prone to the same inaccuracies due to the basis set).
The data in Figure S16 show to relevant aspects: First, the unoccupied states hybridize and broaden much more strongly than the occupied ones (as can be understood from the particularly delocalized character of the LUMO of methane shown in the top panel of Figure S16). Secondly, in the final-state calculations, in which every methane molecule is excited, the DOS is shifted essentially rigidly to higher binding energies by ~4.6 eV. As a consequence, the unoccupied DOS is shifted to the Fermi level and pinning occurs.
Therefore, increasing the charge-transfer distance hardly changes the shift, as is shown for the unoccupied DOS of the 2×2 unit cell in Figure S17. Notably, also the shape of the unoccupied DOS hardly changes upon moving the molecules further away from the surface, indicating that the broadening of the unoccupied DOS in the 2×2 unit cell is really a consequence of inter-molecular (rather than molecule-substrate) interactions.
S27 Figure S16. Occupied Figure S18 for the 4x4 supercell. This effectively reduces the coupling between neighboring molecules and, consequently, the DOS feature that is shifted towards EF for all larger unit cells is a rather sharp peak, as shown in Figures S17 and S19 (with the effect particularly well visible for the 4x4 cell at larger distances). Gaussian including diffuse basis functions as well as with VASP and a plane-wave basis set. These data suggest that pinning at the LUMO would occur at somewhat higher fields in final-state calculations using FHI-aims with a tight basis set compared to equivalent simulations using a plane-wave basis in VASP. The effect is, however, not particularly strong and by no means changes the qualitative picture.

Electrical field for a CREST-type compensation charge
In the main manuscript it has been argued that the spurious collective electrostatic shifts of core-level binding energies in final-state calculations could be avoided by a CREST-like approach 25,26 localizing the compensation charges statically in a charged sheet, rather than putting them into the lowest unoccupied states in the metal next to the Fermi level.
In that case, the charged sheet would have to be placed above the adsorbate layer so that the field between the core holes and the compensation charges is primarily found above the adsorbate molecules and in this way cannot shift the positions of the core levels relative to the metal's Fermi level. A similar effect as with a charged sheet can also be achieved for point charges sufficiently far away from the interface (in the following plot at a distance of 30 Å).
As shown in Figure S20 for two differently sized unit cells, the field of that array of point charges (formed by the core holes and their compensation charges) decays rapidly on the substrate-side of the core holes (i.e. at negative distances), while it quickly approaches a constant value on the other side. This shows that in this way, spurious collective electrostatic shifts of core-level binding energies can indeed be avoided. What the plot, however, also shows is that the decay of the field on the substrate side is much faster than it would be for an isolated point charge. The field of the latter is, however, the field that triggers the screening effects in the metal, that are meant to be described by the final-state calculations. Consequently, a CREST-like compensation charge will result in an incorrect description of the core-hole screening, especially for small unit cells.