Unraveling the Origin of the Repulsive Interaction between Hydrogen Adsorbates on Platinum Single-Crystal Electrodes

Hydrogen adsorption on platinum (Pt) single-crystal electrodes has been studied intensively in both experiments and computations. Yet, the precise origin and nature of the repulsive interactions observed between hydrogen adsorbates (Hads) have remained elusive. Here, we use first-principles density functional theory calculations to investigate in detail the interactions between Hads on Pt(111), Pt(100), and Pt(110) surfaces. The repulsive interaction between Hads on Pt(111) is deconvoluted into three different physical contributions, namely, (i) electrostatic interactions, (ii) surface distortion effect, and (iii) surface coordination effect. The long-range electrostatic interaction, which is generally considered the most important source of repulsive interactions in surface adsorption, was found to contribute less than 30% of the overall repulsive interaction. The remaining >70% arises from the other two contributions, underscoring the critical influence of surface-mediated interactions on the adsorption process. Surface distortion and coordination effects are found to strongly depend on the coverage and adsorption geometry: the effect of surface distortion dominates when adsorbates reside two or more Pt atoms apart; the effect of surface coordination dominates if hydrogen is adsorbed on neighboring adsorption sites. The above effects are considerably less pronounced on Pt(100) and Pt(110), therefore resulting in weaker interactions between Hads on these two surfaces. Overall, the study highlights the relevance of surface-mediated effects on adsorbate–adsorbate interactions, such as the often-overlooked surface distortion. The effect of these interactions on the hotly debated adsorption site for the adsorbed hydrogen intermediate in the hydrogen evolution reaction is also discussed.


Difference between ∆𝑮 diff
ad and ∆ ave ad and their relation with The differential and average adsorption energy (∆ diff ad and ∆ ave ad ) given in Eqs.2-3 in the main text and repeated here, are defined here as in the text by Aravind and Janik, 1 ∆ diff ad =  H * - (-1)H * - ∆ ave ad = where all variables are defined in the main text.
In the Frumkin isotherm, the (potential-dependent) Gibbs free energy of adsorption is assumed to increase linearly with coverage according to 2 Δ ad = Δ non-interact +    ⋅  ⋅ , where  is the interaction parameter,   the Boltzmann constant and  the temperature.In this notation,  is positive for repulsive interaction.Finally, Δ ad and Δ non-interact are the (standard) free energy of adsorption for interacting and noninteracting (Langmuir-type) adsorbates.From the above expression, it becomes clear that  is related to the (constant) slope  of the adsorption Gibbs free energy Δ ad via where  is the slope of Δ ad as a function of coverage.
Considering that we defined two different adsorption energies (∆ diff ad and ∆ ave ad ) in equations (S1a) and (S1b), it may, however, not be immediately clear which of these quantities should enter equation (S2) in order to determine the interaction parameter  computationally.In the following, we derive the Frumkin isotherm to clarify the (approximate) connection between ∆ diff ad and .This connection will become exact in the limit of a large number of possible adsorption sites   and a large number of adsorbates  considered.
We start by considering a system with n hydrogen adsorbates and (m-n) protons in solution.The total Gibbs free energy (including the configuration entropy) of this system is given by The transition from eq.S3a to S3b is performed in accordance with the computational hydrogen electrode (CHE) method. 3For large enough number of adsorption sites   , the configuration entropy can be approximated as 4 where  site is the number of available adsorption sites and  = / site is the coverage.
When plotting  tot () as a function of potential  RHE , in a finite supercell, one obtains plot as shown below.
The transition from a state with  -1 adsorbates to the state with  adsorbates, will occur at the potential  We can now try to simplify Eq. (S5).First, we recognize the left hand side of Eq. (S5) as the expression for ∆ diff ad .In the limit of large numbers, the right hand side can be further approximated as Combing eq (S5)-(S6), we arrive at Comparing this with the Frumkin isotherm within the Nernstian approximation, 5,6 ln Where  is a constant that depends on the equilibrium potential at zero coverage, we find Hence, if the system under investigation can be described via a Frumkin isotherm, Δ  diff ad should be (at least approximately) a linear function of the coverage , and  can be extracted from a linear fit Δ diff ad =  +  between ∆ diff ad and  (at or any other potential) as which is equivalent to Eq. 7 in the main text and  is the slope of ∆ diff ad as a S5 function of .
For H adsorption on Pt(111), we indeed find a nearly perfect linear slope (see Figure . 1 in the main text), justifying the extraction of an interaction parameter  from the slope of ∆ diff ad .

Test calculations for the choice of DFT functional
The average interaction energies at the coverage of 1 ML, ∆ inter ( = 1), were evaluated to check the influence of the choice of density functional on the results.As ∆ inter ( = 1) captures the change in adsorption energy with coverage, its stability (with respect to different functionals used and different parameters used in the DFT setup) is expected to be a good indicator for the convergence of the adsorbateadsorbate interaction (as a function of coverage) that we are interested in in this paper.
As shown in Table S1, the adsorption energy obtained with PBE in the low coverage limit is in better agreement with experimental data.Values obtained for ∆ inter ( = 1) are similar for all three functionals, suggesting that results obtained for the adsorbate-adsorbate interaction (i.e., the change in adsorption energy with coverage) should be robust towards the choice of functional.As can be seen in Figure S1, further converging the k-point grid, the plane-wave energy cutoff or the smearing width, will not change ∆ inter ( = 1) by more than 5% of its absolute value, suggesting that the setting used are sufficiently converged to allow for a meaningful interpretation.

Water bilayer
To investigate the effect of water on H adsorption, we first consider a four-layer slab with a 3×3 hexagonal Pt unit cell and a water bilayer, as shown in Figure S2a-d.The 3×3 supercell was chosen as it is commensurate with a hexagonal water bilayer as considered here.The differential adsorption energy of H adsorption with water bilayer is calculated by where  H * +H 2 O is the Gibbs free energy of the surface with the water bilayer and  hydrogen atoms adsorbed.Other quantities are defined in the main text.
As shown in Figure S3a, the water bilayer has negligible effect on the overall slope of ∆ diff ad as a function of coverage.The water bilayer does, however, influence ∆ diff ad at intermediate coverages, with small energy jumps at 2/9, 5/9 ML and 6/9 ML.We associate these jumps with an energy cost to push the water molecules away from the surface as we see the average Pt -water distance increasing with coverage (see Figure S2c-d).Clearly, these jumps would influence some of the conclusions made in our paper.However, these jumps may easily be overly pronounced in a (static) water bilayer description of the water; in more realistic representations of the adjacent water, this effect should be attenuated.As the average behavior follows the trends observed without a water bilayer, we conclude that the water bilayer has negligible effect on the overall interaction between H ad .The slight increase in differential adsorption energy in the region between 0 ML and 0.6 ML, may, however, explain the slight deviations between the interaction parameter  extracted from cyclic voltammograms in the region above 0 VRHE (corresponding to 0 ML to ~2/3 ML coverage) and the  values extracted from the computational results.

Fully explicit water
We also investigate the water effect on H adsorption in a 4×5×4 orthogonal Pt unit cell with a fully explicit water model using ab initio molecular dynamics (AIMD), as shown in Figure S2e-f.The AIMD calculations were performed using CP2K/ Quickstep package using PBE-D3 as exchange-correlation function, the DZV-MOLOPT-SR-GTH basis set, a plane wave cutoff of 300 eV, a time step of 1 fs, and a Nosé-Hoover thermostat.The average adsorption energy of H adsorption with fully explicit water is calculated by where 〈 Pt+H * +H 2 O 〉  is the time averaged Gibbs free energy of the surface with fully explicit water and  hydrogen atoms adsorbed from our molecular dynamics runs.Other quantities are defined in the main text.The average adsorption energies are calculated and compared here because the AIMD calculations were only performed for several representative coverages (0, 0.2, 0.4, 0.6, 1 ML).
As shown in Figure S3b, the fully explicit water model does not have an obvious effect on the overall slope of ∆ ave ad as a function of coverage -at least not within the error bars achievable within a 15ps molecular dynamics run, when extracting snapshots every 50 fs for analysis.
Since the fully explicit water model still does not show an obvious effect on the overall trend of coverage dependent adsorption energy, and since the interaction parameter  that we find without including water differs only by 30% from the experimentally determined value, 2 we ignore the effect of water in the following.

Configuration effect for H ad on Pt(111)
To investigate the effect of the exact configuration of the adsorbed hydrogens on the surface, several different hydrogen configurations have been tested.In

Relaxed surface
Figure.S5e shows the resulting average adsorption energies Δ ave ad when relaxing the Pt surface as a function of coverage for the configurations considered.Table S2 summarized the numerical values.

Frozen surface
A similar study is performed while keeping the Pt atoms frozen.As shown in Figure .S6, the exact configuration has a much larger influence on the average adsorption energy Δ ave ad in this case.

Electrostatic interaction between H ad
The adsorption of hydrogen atoms at the Pt surface causes charge transfer, as shown in Figure .S8 for the case of two adsorbed hydrogen atoms.This charge transfer can lead to a dipole moment.For a single H atom at the Pt(111) surface, we find a dipole moment of   = 0.013 eÅ from the dipole correction applied in the DFT calculations.
(Note that the dipole correction was not used throughout the paper, but only for this particular test).
In the following, we show that this dipole moment cannot be responsible for the adsorbate -adsorbate interaction observed at low coverages for a frozen surface.To this end, we approximate the adsorbate -adsorbate interaction by only a dipoledipole interaction.In the far-field approximation, the interaction of two (parallel) dipoles in vacuum is given by where  1 and  2 are the two dipole moments,  0 the permeability of vacuum, and  the distance between two dipoles.In the following we consider the interaction energy of one dipole placed at (,) =   2,0 (where   is the x-dimension of the supercell) with all periodically repeated images of a second dipole at (,) = (0,0).
Setting  1 =  2 =   = 0.013 eÅ and considering the size of the supercell, this gives Δ tot dipole-dipole = 0.04 meV, considerably smaller than the remaining adsorbate -adsorbate interaction that can be found from Figure . 1 in the main text.
Several comments are due at this point: First, as shown in Figure .S8, the charge density difference is quite extended in space.A point dipole and the far field approximation will thus be a poor approximation.Second, the charge density difference will not only give rise to a dipole moment, but also to multipole moments, S18 which could be considerably larger, leading to a larger interaction.The point of the exercise above is solely to prove that a dipole -dipole interaction, as often considered for adsorbate -adsorbate interactions, is not sufficient to explain the long-range adsorbate -adsorbate interaction between H adsorbates observed on frozen Pt surfaces.
where  1 and  2 are the Gibbs free energy for x H ad on fcc sites and (x-1) H ad on fcc sites and 1 H ad on top site, respectively.∆ =  2 - 1 .

(=
-1)→ RHE at which the Gibbs free energy of the two states are in equilibrium, leading to  conf () - conf ( -1).(S5)In a finite system, the coverage will change promptly from  -1 to  as the potential decreases from above  (-1)→ RHE to a value below.If  and   are large (for ever-increasing supercell sizes), the coverage  = /  can be changed smoothly, by changing  RHE .

Figure
Figure S1 shows convergence tests for ∆ inter ( = 1) as a function of the number of k-points used to sample the Brillouin zone, the plane-wave energy cutoff energy and the smearing width .The values used to compute final results are highlighted in red.As can be seen in Figure S1, further converging the k-point grid, the plane-wave

Figure S1 .
Figure S1.Parameters convergence analysis in DFT calculations.(a) k-point grid (where  * ) denotes a k-point grid of  ×  × 1) (b) plane wave cut-off energy, (c) smearing width.Unless otherwise mentioned, we use a 4×4×4 orthogonal Pt cell, a 5×7×1 k-point grid, a plane wave cutoff energy of 450 eV, and a smearing width of 0.05 eV.

Figure
Figure S3.(a) Comparison of the differential adsorption energy as a function of coverage for the system without water and with water bilayer (b) Comparison of the average adsorption energy as a function of coverage for the system without water and with fully explicit water.
Fig. S4, we show the different (symmetry inequivalent) configurations for 2 hydrogens adsorbed on a 4x4 surface unit cell.In Fig. S5, we schematically show different configurations considered at higher coverages, namely placing H adsorbates in rows (Figures.S5a and S5b), in a pulk (Figure.S5c) and as far apart as possible (Figure.S5d).In the "row" configurations, once one row is filled, the next H atoms are placed in a second row, that can either be adjacent to the first (Figure.S5a) or as far from the first line as possible (Figure.S5b).

Figure S4 .
Figure S4.Schematic of 3 different configurations for two H adsorbates on Pt(111) surface.

Figure
Figure S5.(a)-(d) Schematic of 4 different configurations for 8 H adsorbates on Pt(111) surface.(e) Average adsorption energy as a function of coverage (ML) for configurations a-d when relaxing the surface.Dotted lines are linear regressions of the plotted data.

Figure
Figure S6.(a) Average adsorption energy as a function of coverage in configuration (b) forming line, (c) cluster, and (d) placing the H ads as far as possible from each .

Figure S8 .
Figure S8.Lattice displacement of the Pt(111) surface during H adsorption. Grey spheres are the Pt atoms, red sphere is the H ad .The arrays show the direction of displacement of Pt atoms, the arrow length indicates how far the atoms move (see text above).

Figure S9 . 7 .
Figure S9.Top view of charge density differences (Δ =  2H * - * - 2H ) for the surface with 2 H ad ( = 2/16 ML) put far apart with each other.The grey and small red sphere are Pt and H atom, respectively.The yellow and blue sphere represent increasing and decreasing charge density, respectively.Red circle indicates the region in which the charge density difference caused by the adsorption of two different H atoms overlaps.

Figure S11 .
Figure S11.Top and side views of charge density differences between Pt surface with  H ad and bare surface, = (a) 1 and (b) 10.Yellow iso-surfaces denote an increase in charge density; blue iso-surfaces a decrease in charge density.

Figure S12 .
Figure S12.Top view of charge density differences between Pt surface with  H ad and Pt surface with ( -1) H ad ,  = (a) 4, (b) 5, (c) 9. Yellow iso-surfaces denote an increase in charge density; blue iso-surfaces a decrease in charge density.

Figure S13 .
Figure S13.Probability for H adsorption occurring on the top sites.

Table S1 .
Calculated adsorption energy at the coverage of 1/16 ML, ∆ ad ( = 1/16), and the average interaction energies at the coverage of 1 ML, ∆ inter ( = 1), as obtained with three different DFT functionals.The computational setup which has been used for these test calculations correspond to: 4×4×4 Pt slab in an orthogonal cell, 5×7×1 k-points, plane wave cutoff energy of 450 eV.The remaining values are equivalent to those given in the methods section in the main text.

Table S2 :
Configurations effect on ∆G ave ad in different coverage of H ad on Pt(111) surface.