What Is the Real State of Single-Atom Catalysts under Electrochemical Conditions—From Adsorption to Surface Pourbaix Plots?

The interest in single-atom catalysts (SACs) is increasing, as these materials have the ultimate level of catalyst utilization, while novel reactions where SACs are used are constantly being discovered. However, to properly understand SACs and to further improve these materials, it is necessary to consider the nature of active sites under operating conditions. This is particularly important when SACs are used as electrocatalysts due to harsh experimental conditions, including extreme pH values or high anodic and cathodic potential. In this contribution, density functional theory-based thermodynamic modelling is used to address the nature of metal centers in SACs formed by embedding single metal atoms (Ru, Rh, Ir, Ni, Pd, Pt, Cu, Ag, and Au) into graphene monovacancy. Our results suggest that none of these SAC metal centers are clean at any potential or pH in the water thermodynamic stability region. Instead, metal centers are covered with Hads, OHads, or Oads, and in some cases, we observed the restructuring of the metal sites due to oxygen incorporation. Based on these findings, it is suggested that setting up theoretical models for SAC modelling and the interpretation of ex situ characterization results using ultra-high vacuum (UHV) techniques requires special care, as the nature of SAC active sites under operating conditions can significantly diverge from the basic models or the pictures set by the UHV measurements.


Introduction
Single-atom catalysts (SACs) present the ultimate limit of catalyst utilization [1][2][3]. Since virtually every atom possesses catalytic function, even SACs based on Pt-group metals are attractive for practical applications. So far, the use of SACs has been demonstrated for numerous catalytic and electrocatalytic reactions, including energy conversion and storage-related processes such as hydrogen evolution reactions (HER) [4][5][6][7][8][9], oxygen reduction reactions (ORR) [7,[10][11][12], oxygen evolution reactions (OER) [8,13,14], and others. Moreover, SACs can be modeled relatively easily, as the single-atom nature of active sites enables the use of small computational models that can be treated without any difficulties. Hence, a combination of experimental and theoretical methods is frequently used to explain or predict the catalytic activities of SACs or to design novel catalytic systems. As the catalytic component is atomically dispersed and is chemically bonded to the support, in SACs, the support or matrix has an equally important role as the catalytic component. In other words, one single atom at two different supports will never behave the same way, and the behavior compared to a bulk surface will also be different [1][2][3].
Looking at the current research trends, understanding the electrocatalytic properties of different materials relies on the results of the physicochemical characterization of these materials. Many of these characterization techniques operate under ultra-high vacuum (UHV) conditions [15,16], so the state of the catalyst under operating conditions and during the characterization can hardly be the same. Moreover, potential modulations under electrochemical conditions can cause a change in the state of the catalyst compared to under UHV conditions. A well-known example is the case of ORR on platinum surfaces. ORR commences at potentials where the surface is partially covered by OH ads , which acts as a spectator species [17][18][19][20]. Changing the electronic structure of the surface and weakening the OH binding improves the ORR activity [20]. Moreover, the same reaction can switch mechanisms at very high overpotentials from the 4e-to the 2e-mechanism when the surface is covered by underpotential deposited hydrogen [21,22]. These surface processes are governed by potential modulation and cannot be seen using some ex situ surface characterization technique, such as XPS.
However, the state of the electrocatalyst surface can be predicted using the concept of the Pourbaix plot, which connects potential and pH regions in which certain phases of a given metal are thermodynamically stable [23,24]. Such approaches were used previously to understand the state of (electro)catalyst surfaces, particularly in combination with theoretical modeling, enabling the investigation of the thermodynamics of different surface processes [25][26][27].
The concept of Pourbaix plots has not been widely utilized for SACs, although there are some examples of their construction [12,28]. However, their use would be exceptionally helpful for understanding the nature of the active sites in SACs under operating conditions and the proper modelling of SACs using computational approaches of different complexity. The latter is specifically related to the fact that the majority of computational models that have been used so far to address he catalytic activity SACs treat SACs as a perfect (single atom + support) combination and do not consider possible changes of the active site due to the potential or pH changes (which are in catalysis, as a rule, rather extreme). Moreover, the use of Pourbaix plots is widespread in electrochemistry and puts the results of DFT thermodynamic calculations in direct connection with the experimental stability of different phases that are present in an electrochemical cell.
In this work, we investigate model SACs consisting of single metal atoms (Ru, Rh, Ir, Ni, Pd, Pt, Cu, Ag, and Au) that have been embedded into a single-vacancy graphene site. Such models have been present in the literature for a while [29]. The incorporation of 3D transition metals, noble metals, and Zn in graphene's single vacancy was studied in detail in Ref. [30]. The reactivity of graphene with a single vacancy (vG) towards the elements of rows 1-6 of the periodic table of elements, excluding lanthanides, is reported in detail in Ref. [31], and the high thermodynamic stability of such systems is observed. Moreover, such systems have also been implemented experimentally and have shown appreciable electrocatalytic activities [32,33]. We start with pristine models of SACs and consider several surface processes, connecting them into Pourbaix plots for given model SACs at the end. We show that the predicted thermodynamically stable states of model SACs change with electrode potential and pH. In fact, the model SACs are actually never pristine, which is the opposite of usual assumptions in the theoretical models of SACs (re)activity that have been considered so far.

Results
To evaluate the stability of different SACs structures under electrochemical conditions, we considered the reactivity of model SACs (M@vG systems) with H, OH, and O. The purpose of this was to estimate which potential regions metal center dissolution (Equation (1)), hydrogen underpotential deposition (UPD, Equation (2)), and the oxidation of metal centers (Equations (3) and (4)) can take place in. To be specific, the considered redox processes were: Once the total energies of the investigated systems were known, and the adsorption energies of the studied adsorbates were determined, it was possible to evaluate standard potentials (E • (O/R)) and to construct the surface Pourbaix plots for the investigated systems (see Section 4 for more details). For reactions (1)-(4), the Nernst equations (at 298 K) were given as:

M@v-Graphene-Formation of SACs
First, we investigated the embedding of Ni, Cu, and Ag and the noble metals Ru, Rh, Pd, Ir, Pt, and Au into the single vacancy site in graphene, i.e., the formation of SACs. When the chosen metal atoms were incorporated into the vacancy of vG, the so-obtained C 3 M sites in the M@vG structures were qualitatively very similar, showing C 3v symmetry in most cases ( Figure 1). In all of the cases, the metal atom protruded from the graphene basal plane, and to a lesser extent, its first three C-neighbours protruded into the plane as well ( Figure 1 and Table 1). The exception to the perfect C 3v symmetry of C 3 M can be found in Ag@vG and Au@vG. Not all M−C bonds have the same length in these systems due to Jahn-Teller distortion (up to the second decimal in Å, Figure 1 and Table 1). From the investigated metals, Ag shows the weakest binding, and Ir shows the strongest ( Table 1). The calculated energies caused by embedding M into the vacancy of vG are in good agreement with available literature reports (Table 1). For the metals belonging to groups 8 and 10 of PTE, we found the total magnetization of M@vG to be equal to zero, while for M from group 11, the total magnetization of M@vG was around 1 µ B (Table 1). Bader charge analysis reveals that some charge is transferred from M to graphene in all the cases (Table 1). While a nearly linear relationship between E emb (M) and the charge transferred from M to graphene was found for Ir, Ru, Ni, Pd, and Au; other investigated elements (Cu, Ag, Rh, and Pt) do not follow this trend. The strongest M binding (Ir) case corresponded to the maximum charge transfer from M to graphene (Table 1).  By comparing the metal embedding energies and the corresponding cohesive energies (experimental data [35], Figure 2), it can be concluded that the majority of the studied metals were less susceptible to dissolution when embedded into vG than the corresponding bulk phase, which is in agreement with our previous findings [36]. The exceptions are Ag and Au, which have lower embedding energies than the cohesive energies of bulk phase (absolute values). Before proceeding further, we note that for the electrochemical applications of SACs, their conductivity must be high. Otherwise, Ohmic losses would affect the energy efficiency of an electrocatalytic process. For this purpose, we investigated the densities of states (DOS, Figure 3) of the studied model SACs. None of the systems show a bandgap, suggesting that all of the studied SACs exhibit metallic behavior. Figure 3. Densities of states for the investigated M@vG systems. Total DOS, carbon, and metal states are given. Plots were generated using the SUMO Python toolkit for VASP [37], and the energy scale is referred to the Fermi level.

H Adsorption (H−M@vG)
The first adsorbate we investigated was atomic hydrogen to explore the possible hydrogen UPD at model SACs. Namely, the bulk surfaces of some of the studied metals show H UPD, such as Pt, Pd, Ir, Rh [38][39][40], as a consequence of the exergonic H 2 dissociation process on these surfaces. Therefore, it is reasonable to expect that at least some of the corresponding SACs could show similar behavior. On the other hand, some other metals, such as Ni, build hydrides, so it is essential to understand the interaction of SAC metal centers with atomic hydrogen.
The calculated E ads (H) ( Table 2) show a relatively wide range of adsorption energies of atomic H on the metal centers of SACs ( Figure 4). Interestingly, the weakest interaction is seen for Ni (which interacts strongly with H in the bulk phase [41,42]) and the strongest is seen for Au (which in bulk interacts very weakly with H [41]). The magnetic moments of SACs are quenched upon H adsorption, but in the cases of Cu and Ru, the magnetic moments arise upon H ads formation.  It is important to consider the geometries of H ads on model SACs. As shown (Figure 3), H ads is formed directly on the metal center in all cases. Furthermore, the H ads formation is followed by reducing a partial charge of the metal center compared to pristine SACs (Table 2), except for in the cases of Ag and Ir, where the situation is the opposite. Based on the obtained results, we can conclude that if H ads is formed on the metal center, the center itself is covered by H and cannot be considered a bare metal site.

OH Adsorption (OH−M@vG)
The OH adsorption energies, referred to as the isolated OH radical, are generally more negative than E ads (H), suggesting a stronger M−OH bond than the M−H bond (Table 3). In all of the systems where OH is bonded directly to the metal center, except for Pd@vG, the partial charge of the metal is lower than in pristine SACs. However, for Cu@vG, we observed an interesting ground state where OH is not bonded to Cu but is instead dissociated and bonded to the carbon atoms adjacent to the Cu center ( Figure 5). This finding is a strong indication that exposing Cu@vG to oxidizing conditions could cause the corrosion of the carbon lattice instead of the oxidation of the metal center.

O Adsorption (O−M@vG)
The studied model SACs bind to the O atom very strongly (Table 4). However, in comparison to OH adsorption and particularly H adsorption, the situation is much less uniform. Ru, Rh, Ir, and Pt SACs bind O directly at the metal center ( Figure 6). Ni and Pd SACs do not bind to O directly, but they do bind at the C atom adjacent to the metal center ( Figure 6). In these cases, the coordination of Pd and Ni by the surrounding carbon atoms reduces from three (pristine SACs) to two, and the C−M−C bridge is formed. For the coinage metals, the metal center coordination numbers are reduced to one (Figure 6), while oxygen atoms are incorporated into the vacancy, resulting in the formation of a pyran-like ring. For these metals, while the system is overall oxidized, the metal center itself is reduced, increasing its partial charge compared to the corresponding pristine SACs (Table 4). In contrast, the metal centers that directly bind O become oxidized as they lose an appreciable amount of charge (Table 4, Ru, Rh, Ir, Pt).

Extensive Oxidation of M@vG (2O−M@vG)
The results presented p until this point indicate that the metal centers and the surrounding carbon atoms in SACs are sensitive to oxidation. While the oxidation beyond Equation (4) is not considered in the construction of the surface Pourbaix plots (for the reasons explained later on), here, we present the results considering the addition of one more oxygen atom to the O−M@vG systems (Table 5, Figure 7). The scenario considered in this section could be operative upon the exposure of SACs to the O 2 -rich atmosphere. As seen from differential adsorption energies (Table 5), O−M@vG systems are prone to further oxidation and bind to O easily. However, this process has devastating consequences on the structure of SACs (Figure 7). In some cases, M can be completely ejected from the vacancy site, while the carbon lattice accepts oxygen atoms. Thus, considering the results presented here, the reactivity of M centers in SACs can be considered both a blessing and a curse. Namely, besides the desired reaction, M centers also present the sites where corrosion starts and, ultimately, lead to irreversible changes and the loss of activity.

Surface Pourbaix Plots for M@vG Catalysts
Using the results obtained for the M@vG, H−M@vG, HO−M@vG, and O−M@vG systems, the surface Pourbaix plots for the studied model SACs were constructed. The construction of the Pouraix plots was completed in several steps. First, using calculated standard redox potentials for the reactions described by Equations (1)-(4) and the corresponding Nernst equations (Equations (R1)-(R4)), the equilibrium redox potentials were calculated for a pH from 0 to 14. Metal dissolution, Equation (R1), is not pH-dependent, but H ads and OH ads formation are, and the slope of the equilibrium potential versus the pH line is 0.059 mV per pH unit in all the cases. Then, the stable phases are identified following the rule that the most stable oxidized phase has the lowest equilibrium potential, while the most stable reduced phase is the one with the highest equilibrium potential. For example, in the case of Ru@vG at pH = 0, the most stable reduced phase is H ads -Ru@vG up to the potential of 0.17 V vs. a standard hydrogen electrode (Figure 8). Above this potential, bare Ru@vG should be stable. However, the potential for the formation of OH ads -Ru@vG is below the potential of the Ru@vG/H ads -Ru@vG couple. This means that the state of the Ru-center immediately switches to OH ads -Ru@vG. The OH ads -Ru@vG phase is the most stable oxidized phase, as it has the lowest redox potential. In the case of Pd@vG, the potential for the OH ads -Pd@vG/ Pd@vG couple is above the one for the Pd@vG/H ads -Ru@vG couple, and there is a narrow pH-potential region where the Pd site is bare. The obtained plots for all of the studied model SACs are presented in Figure 8. Regarding electrochemical applications, we note that all of the systems considered in Sections 2.2.1-2.2.3 are conductive, as seen from the corresponding DOS plots (no bandgap, Supplementary Materials, Figures S1-S3). . We note that in all the cases, metal sites oxidation (OH ads or O ads formation is favored over metal dissolution, and the region of metal dissolution is indicated for completeness, except for Ag, to avoid confusion. While we are aware that the picture presented by Equations (1)-(4) is relatively simple, it is sufficient to capture the essential message provided by the obtained surface Pourbaix plots. As it can be seen, none of the SACs studied here can be considered as being in the pristine M@vG model under electrochemical conditions. Within the water stability window, the metal sites are either covered with hydrogen, or the SAC is oxidized in some way-either at the metal center or the carbon lattice.

Discussion
Based on the presented results, it can be concluded that the actual state of SACs under electrochemical conditions is more complex than in theoretical models. While this conclusion is not surprising, it also points that it is expected that the data obtained under UHV conditions during the SACs characterization could show a completely different picture than the actual one. Our results suggest that considered model SACs are either covered by H ads or are oxidized in some way in the water stability region (formed OH ads or O ads ). Metal dissolution is also possible for some of them in a certain potential-pH range (Ni, Cu, Ag, and Au), but thermodynamically speaking, these metals are anodically protected. Ru, Rh, and Ir centers are either covered by H ads or OH ads . The situation is similar for Ni and Pt. In an aqueous solution, Ni is practically always in the oxidized form (OH−Ni@vG). Pt@vG is covered by H ads at low potentials, while the state switches to the oxidized form at potentials close to the "double layer" region of bulk Pt surfaces [38,43]. It is the same for Pd, where there is a narrow potential region in which the Pd center is bare. However, this region is below the potential of H 2 evolution. Otherwise, Pd is covered by OH ads . In the coinage metal group, Cu behavior is similar to that of Ni. The low stability of Ag and Au at the graphene vacancy ( Figure 2) makes these SACs impractical, and they will not be discussed any further. However, their state switches from H ads -covered to O ads -covered.
The obtained results have implications in terms of both the theoretical modelling of SACs and their practical applications. While the current study is restricted to the metal embedded in the graphene single vacancy site, it should be noted that recently, it was shown that oxygen reconstitutes the Cu-N 2 C 2 site in Cu-based SACs and that Cu-N 2 C 2 -O site is the actual active species of alkaline ORR, while the oxygen reconstitution process is not operative at a low pH [12]. Such oxygen reconstitution is seen here as well ( Figures 5 and 6) for Cu@vG. Moreover, the online ICP-MS showed that the dissolution behavior of as-synthesized Pt SACs with S-containing carbon, as the support is significantly different from that of metallic Pt/C and that the SACs are more stable [16]. However, the S ligands, which stabilize Pt, are prone to oxidation at high potentials (1.5 V vs. RHE), leading to the loss of stability upon the leaching of S and its oxidation. However, the most straightforward confirmation of our conclusions regarding the importance of the realistic nature of SAC probably relates to the in operando characterization of FeN 4 -based SACs under ORR conditions [44]. While this class of SACs is different from the one we investigated, it is of the utmost importance to emphasize that the ORR activities of FeN 4 -based SACs were dictated by the dynamic structure associated with the Fe( 2+ / 3+ ) redox transition and not the static structure of the bare sites. In this specific work, the Fe centers were found to be covered by O ads or OH ads , depending on the potential in the ORR overpotential range [44].
Overall, the results presented here and in previously published experimental findings [12,16,44] indicate that theoretical models for SACs should be carefully set to match realistic operating conditions. Moreover, UHV characterization results must be taken with care when interpreting the electrochemical performance of SACs. In fact, the in situ characterization techniques might be a better choice for extracting SAC properties. Finally, this approach for modelling and evaluating SACs could help us better understand the nature of active sites in these advanced catalysts and point to new strategies for designing single-atom catalysts. To emphasize the importance of considering the state of SACs under realistic conditions, we point to Figure 9. For both Ni@vG and Cu@vG, there is a tremendous impact on the electronic structure upon the oxidation of the metal center according to reaction (1). The well-known relationship between the electronic structure and catalytic activity [45] makes knowing the exact state of the metal centers under electrochemical conditions essential.

Materials and Methods
Graphene with a single vacancy (vG) was obtained by removing one C atom from the graphene plane modelled using a 4 × 4 cell (C 32 ) and relaxing the structure. The 4 × 4 cell was previously confirmed as being large enough to provide valid results for the purposes of this study [30] and Ref. [31]. M@vG systems (C 31 M) were obtained by embedding metal atoms into the single-vacancy site of vG.
The first-principle DFT calculations were performed using the Vienna ab initio simulation code (VASP) [46][47][48]. The generalized gradient approximation (GGA) in the parametrization by Perdew, Burk, and Ernzerhof [49] combined with the projector augmented wave (PAW) method was used [50]. The cut-off energy of 600 eV and Gaussian smearing with a width of σ = 0.025 eV for the occupation of the electronic levels were used. A Monkhorst-Pack Γ-centered 10 × 10 × 1 k-point mesh was used. Chosen metal atoms were placed at the SV site, and during structural optimization, the relaxation of all of the atoms in the simulation cell was allowed. The relaxation proceeded until the Hellmann-Feynman forces acting on all the atoms became smaller than 10 −2 eV Å −1 . Spin polarization was included in all calculations. To include dispersion interactions, which are not accounted for in PBE, we used the DFT+D3 approach of Grimme [51], which corrects the total energy by a pairwise term.
The energy of the metal atom embedding into the single vacancy site of graphene is quantified as its embedding energy:  Vibrational analysis was used to confirm that the relaxed systems are in their stable ground states and to evaluate zero-point energies (ZPE) and vibrational contributions to the entropy (TS vib ). Then, the standard potentials (vs. standard hydrogen electrode) were calculated for the reactions given by Equations (2)-(4) using the total energies of the individual systems (E tot ), zero-point energies, and vibrational entropy contributions (at 298 K), using the computational hydrogen electrode approach [20]. In other words, the equilibrium of the hydrogen electrode is considered: Matching the electrochemical potential of (H + + e − ) to that of H 2 at pH = 0. Further, for each of the competing phases I, the chemical potential (µ i ) was calculated as: The effects of electric field were disregarded, as explained in ref. [25]. To obtain the chemical potential of liquid water, we calculated its chemical potential at the gas phase at 298 K and 1 atm and then corrected it by means of the Gibbs free energy change (∆G) for evaporation under specified conditions. When the chemical potentials for all phases were obtained, the equilibrium potentials for reactions (1)-(4) were calculated by taking the given reaction (Equations (1)-(4)) as a cathode in a hypothetical galvanic cell with a hydrogen electrode as anode, similar to the approach used in [26]. First, the Gibbs free energy changes (∆G) were calculated for each reaction as: The electromotive force of a hypothetical galvanic cell is obtained by dividing calculated ∆G (in eV) with the number of electrons exchanged in the reaction. As the anode is a standard hydrogen electrode, Equation (7), and its standard potential is 0 V, the obtained values are numerically equal to the standard electrode potentials for reaction (2)-(4). For the reaction given by Equation (1), the dissolution potential of M from vG was calculated using the approach described in Ref. [31], considering a hypothetical galvanic cell where one electrode is a massive piece of metal M, while the other is M@vG electrode. For a full description, the reader is referred to ref. [31]. While constructing the surface Pourbaix plots, we kept the concentration of [M z+ ] = 1 × 10 −8 mol dm −3 .

Conclusions
Developing novel electrocatalysts requires an understanding of the nature of the active sites under operating conditions. Our DFT-based modelling of SACs comprised of single metal atoms embedded in graphene monovacancy suggests that under realistic conditions, the nature of metal sites changes, and metal centers are covered with adsorbed H, OH, or O, depending on the electrode potential and pH. We suggest that these changes of the active site architecture under operating conditions must be considered when modelling SACs and when interpreting the results of the ex situ physical and chemical characterization of SACs but should also be taken into account when designing novel SACs.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.339 0/catal11101207/s1, Figure S1: Densities of states for the investigated MG systems with adsorbed hydrogen, Figure S2: Densities of states for the investigated MG systems with adsorbed OH, Figure S3