Surprising Lateral Interactions between Negatively Charged Adatoms on Metal Surfaces

Even something as conceptually simple as adsorption of electronegative adatoms on metal surfaces, where repulsive lateral interactions are expected for obvious reasons, can lead to unanticipated behavior. In this context, we explain the origin of surprising lateral interactions between electronegative adatoms observed on some metal surfaces by means of density functional theory calculations of four electronegative atoms (N, O, F, Cl) on 70 surfaces of 44 pristine metals. Four different scenarios for lateral interactions are identified, some of them being unexpected: (i) They are repulsive, which is the typical case and occurs on almost all transition metals. (ii, iii) They are atypical, being either attractive or negligible, which occurs on p-block metals and Mg. (iv) Surface restructuring stabilizes the low-coverage configuration, preventing atypical lateral interactions. The last case occurs predominantly on s-block metals.


S5 Adsorption Induced Work Function
Changes S7 S1 Description of the Simple Ionic Model The "simple ionic model" was derived in our previous publication 1 and here we briefly explain its essence. The model is based on the electronic structure analysis of the O/Al(111) system, which reveals that O-Al bonding is ionic and furthermore that the excess electron charge on adatoms mainly comes from the nearest neighbor metal atoms. 1 The last observation is exploited in the simple ionic model, where the adatom gets the electron charge exclusively from the nearest neighbor metal atoms, such that each neighboring metal atom contributes proportionally, as schematically shown in Figure S1. The simple ionic model therefore consists of an ionic bilayer of adatom-anions/metal-cations and can be described by N ions in the unit-cell at positions {τ i } N i=1 and charges {q i } N i=1 . The interaction energy is then obtained by summing the pairwise Coluomb interactions among the ions in the infinite adatom/metal bilayer, i.e.: where {R} are the lattice vectors of a two-dimensional lattice. The role of the (1 − δ i,j δ R,0 ) term is to omit the interaction of an ion with itself (i = j and R = 0, where R = |R|). The infinite lattice sum in two-dimensions, ∞ R =0 (· · · ), can be evaluated by explicitly calculating it within the cutoff radius R cut , whereas beyond R cut it is approximated by an integral, in particular:  Figure S1. Schematic of the electron charge distribution among adatoms and nearest neighbor metal atoms in the simple ionic model. (a) The charge on adatoms is set to −q ad (labeled as "−"), whereas the counter-charge is distributed proportionally to the nearest neighbor metal atoms, i.e., for an adatom with n metal neighbors each of them donates q ad /n electrons to the adatom. If a metal atom has no adatom neighbors then it remains charge-neutral, but if it has m adatom neighbors then it donates q ad /n electrons to each, thus in total mq ad /n to all of them. (b) Charge distribution for (2 × 2) and (1 × 1) adatom overlayers on a square lattice of metal atoms, which is compatible with bcc(100), fcc(100), and sc(100). Unit-cells are indicated in red.
where A is the area of the unit-cell and z i is thê z-coordinate of the atomic position τ i , i.e., τ i = (x i , y i , z i ). Note that the larger is the R cut , the smaller is the last sum (− 2π A N i,j · · · ), which is due to charge neutrality and consequently the cancellation between its terms. In particular, the last sum scales as (2π/A)µ 2 z R −1 rcut , where µ z is theẑ-component of the dipole of the ions in the unit-cell, µ z = N i=1 z i q i . Hence:  Figure S2. The prediction of the simple ionic model for the (1 × 1) and (2 × 2) adatom overlayers above a square lattice of metal atoms, which is compatible with the bcc(100), fcc(100), and sc(100) surfaces. Notice that below the critical height the high-coverage (1×1) overlayer is more stable than the lower-coverage (2 × 2) overlayer. (a) The dependence of the interaction energy and the critical height on the lattice parameter for a = 1.6Å and a = 5.3Å. Notice that the (1 × 1) over (2 × 2) preference below the critical height decreases by increasing the lattice parameter, while concomitantly the critical height becomes larger. (b) The dependence of the interaction energy on the charge of the adatom for q = −0.5 and q = −1.0. Because the energy depends quadratically on the charge, the (1×1) over (2×2) preference below the critical height decreases with a decrease in charge magnitude. Figure S2 depicts the results of the simple ionic model for the (1×1) and (2×2) adatom overlayers over a square lattice of metal atoms, with the interaction S2 energy shown as a function of adatom height above the surface. The figure illustrates the dependence of the interaction energy on the lattice parameter (Figure S2a) and on the adatom charge ( Figure S2b). The most important result of the ionic model is that there exist a critical adatom height below which the highcoverage (1 × 1) overlayer is more stable than the low coverage overlayer (currently only the (2 × 2) overlayer is considered for low coverage, but in our previous work 1 we considered even lower coverages). This effect is referred to as "stabilization" in the following. The critical height obviously depends on the lattice parameter due to geometric reasons, i.e., the larger is the lattice parameter, the larger is the critical height. Concomitantly with the increase of the lattice parameter the extent of stabilization obviously decreases ( Figure S2a). As for the dependence on the adatom charge, the stabilization obviously increases with the magnitude of the adatom charge ( Figure S2b), because the energy depends quadratically on the charge. Some of these dependencies can be easily understood from eq (S1).

S2 Specification of PAW potentials
We used PAW potentials from the pslibrary (https://dalcorso.github.io/pslibrary), mainly from version 1.0.0. For a few metals potentials from older versions of the pslibrary were used, in particular:

S3 Binding and Adsorption (Free) Energies
Binding energies were calculated with the following equation: where X stands for the adatom (either N, O, F, or Cl) and E X /slab , E X , and E slab are DFT total energies of the adatom/slab system, standalone adatom, and bare slab, respectively. The adsorption energy, measured with respect to isolated X 2 molecules, is defined as: where E X 2 is the total energy of the standalone diatomic molecule. E ads can be obtained from E b as: where D X 2 is the bond strength of the X 2 molecule: By plugging the DFT calculated values of D X 2 into eq (S6), the following relations are obtained: Note that the so calculated E ads is normalized to a single adatom. As a first approximation, dissociative chemisorption of X 2 is exothermic when E ads < 0. As to get a better estimate of whether a given adatom overlayer is themodynamically stable with respect to gas-phase X 2 molecules, we also calculated adsorption free energies at a partial pressure of p = 1 atm and T = 298.15 K. These calculations were performed by following the approach described in the Supporting Information of our previous publication. 2 To this end, all degrees of freedom were relaxed and vibrational frequencies were evaluated at the Γ q-point using the PHonon code 3 from the Quantum ESPRESSO distribution. 4 To account for the breakdown of the harmonic oscillator model for low-frequency vibrational modes, we utilized a quasi harmonic approximation of Cramer-Truhlar 5 within which the vibrational frequencies below a threshold value corresponding to hν = 1 2 kT are raised to the value corresponding to 1 2 kT ; for T = 298.15 K we used a threshold value of 100 cm −1 .
For convenience, let us write the free energy as: where E 0 is the Kohn-Sham total energy at zero temperature without the zero-point energy (ZPE). For the gas-phase molecules G corr is given by: where E trv (T ) and S trv (p, T ) are translationalrotational-vibrational thermal energies (with ZPE included) and entropies, respectively. For the oxygen molecule, which is in a triplet ground state, also the electronic degrees of freedom were taken into account, S elec = k ln (3). For slabs only the vibrational contribution to thermal energy and entropy was taken into account, whereas the pV term was neglected, i.e.: where E vib (T ) and S vib (T ) are the vibrational thermal energy (with ZPE included) and entropy, respectively. Configurational entropy of slab systems was not taken into account, because only fully ordered structures are considered (perfect bare surfaces and (1×1) and (2×2) adatom overlayers). The adsorption free energy was then calculated as: where: (S13) ∆G corr values were calculated for several adsorption systems at p = 1 atm and T = 298.15 K (Table S1) and the values span the 0.3 ± 0.1 eV range. Hence, as a first approximation, we can write: at p = 1 atm and T = 298.15 K.
This result can be rationalized as follows: for gasphase molecules the calculated roto-translational contributions to the Gibbs free energy at room temperature and partial pressure of 1 atm are −0.53, −0.54, −0.56, and −0.62 eV for N 2 , O 2 , F 2 , and Cl 2 , hence the corresponding contributions to the adsorption free energy are 0.26, 0.27, 0.28, and 0.31 eV, respectively. These numbers are not far from the ∆G corr values tabulated in Table S1. This implies that vibrational contributions to the adsorption free energy are small (below 0.1 eV) due to cancellation between reactants and products.

S3.1 Comparison of methods: PBE vs. PBE+U+D2
Unless explicitly stated otherwise all the results reported herein refer to the PBE functional. The PBE S1: ∆G corr values, calculated at 1 atm and 298.15 K, for several adsorption systems. For a given adatom stronger bonding to the surface usually implies a larger ∆G corr value, due to ZPE in the adsorbed state (note that among the three presented surfaces, adsorption bonding is the strongest on Al, followed by Na, whereas the bonding on Cu is the weakest).   functional was reported to be the "best compromise" among the family of gradient-corrected functionals for the description of binding energies, surface energies, lattice constants, and adsorption energies for late transition metals; 6 the same study also argued that "for metals and metal surfaces, hybrid functionals are hardly a major step forward". As to illuminate how the choice of the PBE functional may affect the results, we also made a few tests with the PBE+U+D2 method, where +U implies the simplified version of the GGA+U method of Cococcioni and de Gironcoli 7 and +D2 stands for the DFT-D2 semi-empirical van der Waals correction of S4 Grimme. 8 To this end, we performed PBE+U+D2 calculations for the four adatoms on Fe(100) and Fe(110) surfaces. Iron was chosen because it has been considered as a benchmark d metallic system for the LDA+U method in the seminal papers of Anisimov and Gunnarsson 9 and Cococcioni and de Gironcoli. 7 The value U = 2.2 eV was taken for Fe, as suggested by Cococcioni and de Gironcoli. 7 The comparison between PBE and PBE+U+D2 calculated adatom binding energies is provided by Table S2. While the absolute values differ by up to about 0.3 eV (PBE+U+D2 gives weaker bonding for N adatom and stronger bonding for F and Cl adatoms), qualitatively the two methods give compatible results. Most importantly the differences between the (1 × 1) and (2 × 2) binding energies (∆ E b ), which are the essence of the current study, are quite comparable between the two methods; PBE+U+D2 gives slightly smaller ∆ E b values. Because the main purpose of the current study is the identification of attractive interactions, slightly larger PBE ∆ E b differences may suggest that PBE is a bit more conservative and safer to use for the current purpose.

S3.2 Tabulated Binding Energies
The binding energies for adatom overlayers on bcc, fcc, and hcp metal surfaces are tabulated in Tables S3 to S7 and plotted in Figures S7 to S11; these figures also consider if adatom overlayers are stable with respect to gas-phase X 2 molecules (X = N, O, F, or Cl). The E b values of sc (simple-cubic) systems are given in Table S8, but they are not plotted, because only Bi is "considered" to crystallize in this lattice type. Note that Bi and some other investigated metals that crystallize in more "exotic" lattices were instead modeled-to simplify the calculations-by one among the bcc, fcc, hcp, or sc lattice types, as described in the main text.

S3.3 Stability of Adsorption Sites
The adatoms mainly adsorb in hollow sites, however, the number of exceptions is significant. In such cases hollow sites are either less stable than bridge, top, or subsurface interstitial sites or they are even unstable. The most stable identified sites for each case are indicated in Tables S3-S8, which  fcc-and hcp-hollow, whereas for bcc(110) surfaces hollow sites can be either three-fold or fold-fold coordinated (cf. Figure S3).
We begin our analysis by noticing that our results are in line with those reported by Zhu et al. 10 for halogen adatoms (see the comparison in Figure S12). Namely, both F and Cl prefer the top site on Al(111) at the lower coverage. The top site is also the preferred site for (2 × 2) adlayer of F on Ir(111) and for (2×2)-F and (1×1)-Cl on Pt(111). Note that most of these exceptions have been reported by other authors as well. 11 As for bcc metals the top site is preferred for F on W(110). However, we find that F prefers the hollow site on Mo(110) and not the top site as reported by Zhu et al. 10 The site preference for F and Cl on hcp metals is also reproduced. The calculated E b values for (2 × 2) layers of F and Cl are compared to those reported by Zhu et al. in Figure S12. For Cl the average difference between the two sets of E b values is −0.01 ± 0.09 eV, whereas for F the average difference is 0.21 ± 0.13 eV.
In addition to the already documented site anomalies, we found that (2 × 2)-F prefers to adsorb on the top site of Os(001). (1 × 1)-F and (1 × 1)-Cl also adsorb on the top site of Ru(001) and Re(001). The top site is also the most stable for F on W(110) and for F on Bi(100) at low coverage.
Finally, we comment on the three-fold and four-fold hollow site occupation on bcc(110) surfaces. Adatoms adsorb predominantly in the three-fold hollow site; the exceptions are listed in Table S4. An example of three-fold and four-fold hollow sites on bcc(110) is shown in Figure S3. Figure S4.
The unoccupied surface area (A usa ) is calculated by subtracting the area occupied by the metal cation from the area of the (1 × 1) unit-cell. For the radius of the metal cation, R c , the average value of cationic radii for all coordination numbers of the lowest cationic oxidation state is used. Radii were taken from Shannon and Prewitt. 12

S4 Surface Restructuring
In addition to attractive, negligible, and repulsive lateral interactions we also identified another possibility, when the two conditions of the simple ionic model for the attractive lateral interactions are met (ionic adatom-surface bonding and low height of adatoms). In particular, adatom induced surface restructuring stabilizes the (2 × 2) overlayer and makes it more stable than the (1 × 1) overlayer.

S4.1 Unoccupied Surface Area
As an approximate ad hoc criterion for which adatom/metal pairs restructuring can be expected, we defined a quantity termed unoccupied surface area (A usa ), whose calculation is schematically illustrated in Figure S4.

S4.2 Restructuring Quotient
In order to quantify the extent of surface restructuring for each adatom/metal pair, we defined the restructuring quotient (f res ) by eq (2) in the main text. The way the restructuring quotient is calculated is schematically illustrated in Figure S5. It is evident from this figure that f res pertains to adatoms adsorbed in hollow sites, hence it is considered only for (2×2) overlayer structures with adatoms adsorbed in hollow sites. Note, however, that the hollow sites are not always the stablest (cf. Table S3-S8). We defined the surface to be restructured only when f res is significantly below 1; we arbitrarily set f res ≤ 0.9 as the criterion for surface restructuring. restructuring area of (1×1) cell restructured area , whereas A R is the area enclosed by the nuclei of the four metal atoms nearest to the adatom forming a (2×2) overlayer (right panel).
It is evident from this definition that f res refers only to adatoms adsorbed in hollow sites. The restructuring quotient was used to estimate the degree of restructuring, in particular, we used f res ≤ 0.9 as the criterion for surface restructuring.
We should comment on how the restructuring quotient was calculated for close-packed bcc(110), fcc(111), and hcp(001) surfaces. In particular, on bcc(110) the adatoms were found in two distinct hollow sites, i.e., three-fold and four-fold hollow sites (see Figure S3), whereas on fcc(111) and hcp(001) both fcc-and hcp-hollow sites are three-fold coordinated. For all the three-fold hollow sites, f res was estimated by considering the area spanned by the nuclei of the three metal cations nearest to the adatom, whereas for the bcc(110) four-fold hollow site the area spanned by the nuclei of the four nearest metal cations was taken into account.
The obtained values for the restructuring quotient of the (2 × 2) overlayer structures with adatoms adsorbed in hollow sites are plotted on the right-hand side of Figures S7 to S11. The plots clearly show that the extent of restructuring is the greatest for N and O adatoms on bcc(100) surfaces, where alkali metals stand out the most. However, on bcc(100) surfaces of alkali-metals even more stable structures exist that are characterized by the adatoms located in either subsurface octahedral (O and N) or tetrahedral (F and to lesser extent Cl) sites. These two sites are S6 (1×1)-O@Na(100) (2×2)-O@Na(100)  (100). For the (2×2) overlayer the Na atoms closest to the adatom move toward it so that Na 4 O islands form, whereas for the (1 × 1) overlayer this is not possible due to translational symmetry. (b) However, on alkali-metal bcc(100) surfaces even more stable (2 × 2) structures can form by adatoms migrating below the surface bridge sites into subsurface interstitial sites, either "octahedral" (O and N) or "tetrahedral" (F and to lesser extent Cl); note that bcc lattice does not have proper octahedral and tetrahedral interstitials. The difference between the two sites is in the height of the adatom. For octahedral site the adatom is at the height of subsurface layer and octahedron consists of 2 surface, 2 subsurface, and 2 sub-subsurface metal atoms, whereas for tetrahedral site the adatom is in between the surface and subsurface layer and tetrahedron consist of 2 surface and 2 subsurface metal atoms. For (2 × 2) structures, the occupation of interstitial sites, in particular for the octahedral site, is characterized by attractive displacement of neighboring metal atoms toward the adatom. (c) Such attractive pulling of metal atoms as to form octahedron is not possible for high-coverage (1 × 1) overlayers due to translational symmetry, hence the adatom is instead located in a deformed five-fold coordinated site.
shown schematically in Figure S6b. Indeed, on the basis of f res and A usa values, calculated with eqs (2) and (3) in the main text, we can conclude that when the f res quotient (for structure with adatoms in hollow sites) becomes too small or when the A usa factor is too high, the adatoms will migrate into subsurface interstitial sites, alkali-metal bcc(100) surfaces being a notable example.
In the case of fcc metals restructuring occurs for N, O, and F on the Ca(111) and Sr(111).

S5 Adsorption Induced Work Function Changes
In some cases, such as N, O, and F adatoms on Mg(0001) the lateral attractive interaction between adatoms are accompanied by an adsorption induced decrease of the work function. 13 It should be noted, however, that a decrease of the work function is not required for attraction to emerge, as evidenced by Figure S13, which shows the adsorption induced work function change for all the "attraction cases" currently identified. Among 28 such identified cases, work function reduces for only 13 of them. Figure S14 plots the adsorption induced work func-S7 tion changes for all the considered adatom overlayers and Figure S15 shows the experimental work functions for either closed-packed or polycrystalline surfaces of the 44 metals considered in this study.    Figure S7. Left: binding (E b ) and adsorption (E ads , right ordinate) energies of adatoms on bcc(100) surfaces for (1×1) and (2×2) overlayers. The dark-red line indicates the zero adsorption energy, which equals to − 1 2 D X 2 , cf. eqs (S6) and (S8), whereas the brown band indicates the estimated zero value of G ads with the error-bar of ±0.1 eV, cf. eq (S14), at p = 1 atm and T = 298.15 K. Note that for O ads , F ads , and Cl ads the E ads zero lines are not shown, because their E ads ranges are fully exothermic. Right: restructuring quotient (f res ) of bcc(100) surfaces for (2 × 2) adatom overlayers adsorbed at hollow sites. Cases where a subsurface interstitial site is the most stable are indicated by orange dots. The red dash-dotted line at f res = 0.9 indicates the adopted threshold for the restructuring. Lines are drawn to guide the eye. S9   Figure S8. As in Figure S7, but for bcc (110) Figure S9. As in Figure S7, but for fcc(100) surfaces. Note that for F ads the zero E ads line is not shown, because its E ads range is fully exothermic. For the restructuring quotient, a few cases where subsurface interstitial sites are more stable than hollow sites are indicated by orange dots. S11  As in Figure S7, but for fcc(111) surfaces. Note that for F ads the zero E ads line is not shown, because its E ads range is fully exothermic. S12  ( 1 × 1 ) Figure S11. As in Figure S7, but for hcp(001) surfaces. Note that for F ads the zero E ads line is not shown, because its E ads range is fully exothermic. S13       Adsorption induced work function change (∆Φ) of (1 × 1) adatom overlayers for all identified cases of attractive lateral interactions (∆ E b < −0.1 eV), whereas adsorption induced work function changes of all considered adatom overlayers are shown in Figure S14.    Where available the work function for the close-packed surface is plotted, otherwise the polycrystalline value is used. The lines are drawn to guide the eye. S18