INTRODUCTION

As a new frontier catalytic material, supported single-atom catalysts (SACs) have attracted extensive attention in the recent decade since the concept of "single-atom catalysis" was first introduced by Zhang et al. 1 With the high atom utilization and remarkable catalytic activity, SACs exhibit great potential for applications for electrochemical energy conversion and storage processes, such as hydrogen evolution reaction (HER), 2 oxygen evolution reaction (OER), 3 carbon dioxide reduction (CO 2 RR), 4,5 nitrogen reduction reaction (NRR), 6,7 oxygen reduction reaction (ORR), [8][9][10] etc. Despite numerous studies focusing on the design, synthesis, and application of SACs, their commercialization is still limited due to the pursuit of better activity and selectivity at the expense of stability, which is a more critical factor hindering large-scale applications. [11][12][13] The instability of SACs is attributed to the high surface energy of metal single atoms prompting them to aggregate into more stable nanoparticles during conventional preparations, 14 such as the coexistence of Pt single atoms and nanoparticles, which often occurs in conventional substrates of Al 2 O 3 , SiO 2 , and TiO 2 . 15 Furthermore, SACs can undergo structural transformations into clusters under working conditions. 16 For instance, Fredeŕic et al. found that FeN 4 C 12 and FeN 4 C 10 moieties possess different stabilities with increasing operation time. The number of active sites for the former is reduced due to the formation of ferric oxides, resulting in the deactivation of the catalyst, while the number of catalytic active sites for the latter remains unchanged. 17 In addition, a more interesting structural transformation occurred on Cu SACs during CO 2 reduction. Fontecave et al. 18,19 reported that single-atom Cu 2+ dispersed in the N-doped carbon material, under a certain range of negative electrode potential, could be transformed into Cu 0 small clusters. In turn, it can be restored to Cu 2+ in the initial atomic dispersion state after releasing the applied potential or switching to a positive potential. Similar dynamic behavior was also observed by Xu et al. 20 It is even clearer that copper-based metal-organic frameworks with well-defined structures showed reversible dynamic reconstruction during working conditions. These transformation behaviors not only do not deactivate the catalysts but also act as true catalytic active sites, facilitating the efficient reduction of CO 2 to ethanol or methane. 21 However, the current evaluation of the stability of SACs is still mainly focused on thermodynamic "static" calculations, including formation energy, 22 cohesive energy, 23 binding energy, 2 diffusion barrier, 24 etc. Although molecular dynamics simulations can also be used to verify the stability of SACs, the most relevant studies do not consider the effects of kinetics, pH, working potential, and concentration of reactants or key intermediates on the structural stability in complex and realistic catalytic reactions. [25][26][27] Therefore, exploring the "dynamic" evolution of the structural and electrical properties of SACs under working conditions is a great challenge and is essential.

Herein, we construct a solid-liquid interface, namely, copper single-atom nitrogen-doped graphene (Cu-N 4 -C) with explicit water layers, to study the dynamic stability of the Cu single-atom catalyst under working conditions using our recently developed constant-potential hybrid-solvation dynamical model (CP-HS-DM). 28 The calculation results demonstrate that the more negative the electrode potential, the stronger the adsorption of H at the N site, thus inhibiting the production of H 2 , which is consistent with the experimental result. 18 Moreover, it is found that the adsorption of H is a vital driving force for the leaching of the Cu single atom (SA) from the catalyst surface during electrolysis, which can reduce the kinetic barrier to Cu SA leaching compared to the absence of H adsorption. After a portion of Cu atoms leach from the surface, they form transient steady-state Cu small cluster structures with adjacent and tethered Cu atoms, thus facilitating the reduction of CO 2 to ethanol. When the reaction is over, or the applied potential is adjusted to a positive value, the oxidation of Cu is greatly accelerated by OH • and then redeposited, and finally, the Cu clusters are restored to the initial atomically dispersed state, completing the reconstruction cycle of the copper catalyst. Our work provides a new fundamental understanding of the dynamic stability of SACs under working conditions.

COMPUTATION METHODS

Structure optimization and ab initio molecular dynamics (AIMD) were carried out using density functional theory (DFT) as implemented in the Vienna ab initio Simulation Package (VASP 5.4.4). 29,30 The electronic exchange and correlation interactions were described by the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA). 31,32 The projector augmented wave (PAW) method 33 was applied to account for the core-valence interactions and the cutoff energy of the plane-wave basis was set to be 400 eV in the relaxation, while 300 eV was used in the AIMD simulations at 300 K. The convergence criteria of energy and force are 0.01 eV/Å and 10 -5 eV for relaxing structures, respectively. The PBE-D3 method 34 was employed to correct van der Waals interaction of water-water and water-substrate. A 3 × 3 × 1 gamma-centered k-mesh was used in the relaxation, while the gamma point of the Brillouin zone with no consideration of symmetry was adopted in the AIMD simulations. The time step in AIMD was set to be 0.5 fs. The dynamic stability of the catalyst structures under operating potential was evaluated by CP-HS-DM. 28 The net electronic charges were balanced by the ionic charges in the implicit solution, and thus, the total system remains charge neutral. The compensation charges were added as point charges, following the Poisson-Boltzmann distribution. This has been implemented to

Journal of the American Chemical Society

VASP as VASPsol. 35 More details can be found in the Supporting Information.

RESULTS AND DISCUSSION

The prerequisite of Cu SA leaching from the surface to form a small Cu cluster during CO 2 RR is that the Cu-N bond becomes weak and easy to break. Therefore, the possibility of Cu SA leaching is evaluated in thermodynamics and kinetics separately. The free energies of the leaching of Cu single atoms anchored on the substrate to form solvent cations Cu 2+ (aq) are calculated for determining the thermodynamic stability of the Cu single-atom catalyst. Similar to the computational hydrogen electrode method, 36 the free energy of the Cu ion G(Cu 2+ (aq)) is obtained from the experimental standard hydrogen electrode U 0 and the calculated free energy of the bulk metal G(Cu(s)) as follows

Cu (aq) 2e

Cu(s)

U e e SHE 0

where μ e is the electron energy, μ SHE is benchmarked to be -4.6 eV for VASPsol, 35 and U 0 = 0.34 V. The leaching process is as follows

Thus, the corresponding reaction free energies at a constant potential U can be calculated as Figure 1a shows the free energies of Cu SA leaching from the surface to form Cu 2+ (aq) at different electrode potentials from 0 to -1.5 V vs reversible hydrogen electrode (RHE). Although the free energy decreases to 2.09 eV at -1.5 V vs RHE, it is still thermodynamically difficult to leach from the surface. We further evaluate the kinetic possibility of Cu leaching from the surface by CP-HS-DM. Here, the referenced electrode potential is -1.2 V vs RHE, since the highest faradaic yield of ethanol can be obtained at this corresponding potential. 18 A solid-liquid interface model is constructed in Figure S1. In this structure, an inert He atom is introduced as a rivet and fixed at a distance of 18.28 Å directly above the Cu atom, using it to calculate the kinetic activation energy of Cu SA leaching from the N 4 -C site. The chosen reaction coordinate in the "slow-growth" approach is that the distance between Cu and He atoms is shortened (see Figure S3a). For more visualization, the horizontal coordinates in Figure 1b are converted to the distance between dynamic Cu and initial Cu. The free energy continues to increase as Cu moves away from the surface, reaching up to 4.74 eV at the end of the reaction. This result indicates that the four Cu-N bonds are difficult to break, and Cu SA cannot be leached from the surface at room temperature. The dynamic evolution of the bond length of Cu-N with the Cu SA leaching process is tracked in Figure 1c. Note that one water molecule close to the N (3) atom undergoes hydrolysis when the Cu-N (3) bond is broken and the corresponding configuration at the end of the reaction is that the Cu atom is only coordinated with one N atom and one -OH group, which is in accordance with the experimental EXAFS spectra data. 18 Figure 1b shows the change in the number of electrons with the evolution of the structure. Compared with the initial structure, the net charge number of the final structure increases to about 0.7 e -due to the formation of OH -. Bader charge analysis of the representative snapshots in Figure 1b is shown in Figure 1d, including Cu SA, substrate, water, and He atom. Obviously, the He atom does not gain or lose electrons in the whole process, which has a negligible effect on the charge of the whole system.

Interestingly, the number of electrons lost by Cu decreases from 0.9 e -in the initial structure to 0.56 e -in the final structure, indicating that Cu SA leaching from the surface tends to be metallic Cu, which is in agreement with the experimental observation that the Cu oxidation state changes from +2 to 0. 18 Therefore, both thermodynamic and kinetic results indicate that it is almost impossible for the Cu SA to leach directly from the surface.

However, the existence of a structural transformation between atomically dispersed Cu +2 and metallic Cu 0 small clusters through X-ray absorption spectroscopy (XAS) characterization in operando CO 2 electrolysis conditions has indeed been confirmed. 18,20,21 What is the driving force for the leaching of Cu SA from the N 4 -C site with the strong chelating capacity to form Cu clusters? Does H 2 as a competitive product for CO 2 reduction affect structural transformation? To solve these confusions, the H adsorption behavior under the CO 2 reduction potential is explored. Figure 2a,b shows free energies as a function of the applied electrode potential and configurations of xH adsorption (x = 1, 2, 3, and 4), respectively. As the electrode potential changes from zero to negative, the corresponding free energies of H adsorption convert from positive to a negative value, indicating that the adsorption of H at the N sites is thermodynamically converted from unfavorable to favorable. This transformation is ascribed to the fact that the potential becomes more negative, resulting in more electron accumulation on the catalyst surface, thus promoting the adsorption of H + (see Table S1). Note that the calculations of free energy based on the conventional firstprinciples method show that the H adsorption on N (x) (x = 1, 2, 3, and 4) sites is thermodynamically unfavorable, as shown in Figure S6. These results adequately demonstrate that it is critical to consider the surface charge effect in electrochemical reactions. 38,39 More interestingly, our calculation gives a ΔG H value close to zero at U RHE = -1.0 V, which explains why the faradaic yield of H 2 as the major reaction product is the highest at -1.0 V in the experiment. 18 At U RHE = -1.2 V, the hydrogen evolution reaction (HER) is effectively suppressed due to the significant enhancement of H adsorption. We speculate that the adsorption of H may be an important factor affecting Cu desorption at a negative potential (which will be discussed later). The source of H is further discussed. There is only one possibility, that is, hydrolysis. Therefore, the kinetic barrier for the hydrolysis of the first H 2 O to form H* and OH -is calculated (see Figures 2c andS3b). The criterion for the end of this reaction is that the difference between the N-H and O-H bond lengths is 1.00 Å in Figure S7a, and the energy barrier is 0.64 eV at -1.2 V, indicating that this reaction easily occurs at room temperature. Note that the difference in the number of electrons before and after the reaction is about 0.8 e -, which confirms the production of OH -(see Figure S7b).

To have an idea of how much the randomness of the AIMD can impact the results, incremental speeds of 0.0004 and 0.0012 are considered for the hydrolysis of the first H 2 O, as shown in Figure S8. The calculations show that the energy barrier fluctuates within ±0.05 eV, indicating that the accuracy of the calculations employed in these AIMD simulations is acceptable. Although the energy barrier for H adsorption is relatively low, the energy of the final state is higher than that of the initial state. One may wonder if removing the constraint can automatically reverse the process to the initial state, so an unconstrained AIMD simulation of 2 ps is performed. The calculated results show that it takes 1.37 ps for the adsorbed H to return to the solution (Figure S9a). H returns because there is a OH -in the aqueous solution. Under real electrochemical conditions, OH -may diffuse into the bulk of the solution and reduces the free energy, resulting in adsorbed H not returning to the solution. 40,41 Since the structure of the catalyst contains four N atoms, it can provide four active sites for H adsorption, and the activation energies of the remaining three N sites occupied by H also are calculated, which are 0.62, 0.46, and 0.67 eV (Figures 2d andS10). The average barrier for the hydrolysis of H 2 O to form *H and OH -is about 0.60 eV. Therefore, we verify from both thermodynamics and kinetics that H from the H 2 O molecule can be adsorbed to the N sites at

Next, we try to answer the remaining question above: how does the adsorption of H affect the transformation of the catalyst structure during CO 2 reduction? The free energy (ΔG) of Cu 2+ (aq) formation during Cu SA leaching under 1H and 2H adsorption is calculated and shown in Figure 3a. With the adsorption of H, the binding strength between Cu SA and the substrate is weakened, leading to a thermodynamically unfavorable to favorable leaching process of Cu SA. Back to the dynamics, the kinetic barrier for the Cu SA leaching process from the Cu-N 4 -C substrate is re-evaluated, where one N site is adsorbed by one H. An activation energy of 0.70 eV is taken from the energy of the 2500th step (1.25 ps) in Figure 3b, corresponding to the coordination of Cu with one N atom and one H 2 O molecule. This is consistent with the experimental observation that Cu coordinates with light atoms (C, N, O) during CO 2 RR. We further perform augmented sampling of the snapshot structures in Figure 3b and an AIMD simulation of 2.5 ps for each structure in a fixed collective variable to obtain a kinetic energy barrier of 0.66 eV (see Figure S11), which is consistent with the slow-growth results. In addition, an unconstrained AIMD simulation of 1.5 ps is performed, and it is observed that the Cu single atom is not returned to the initial state, indicating that the final state is a relatively stable state (Figure S9b). Compared with a pure Cu-N 4 -C surface (4.74 eV), the adsorption of H significantly promotes the leaching of Cu SA. As shown in Figure 3c, when the constrained MD reaches 2.25 ps, the Cu atom is completely detached from the surface, dissolved in an aqueous solution and attached to two H 2 O molecules. Unexpectedly, Figure 3d shows the dynamic evolution of the bond length of Cu-N 4(x) (x = 1, 2, 3, and 4), indicating that the Cu SA spontaneously leaches from the surface in a short time (∼300 fs) under the condition of 2H coadsorption. The obtained final configurations are composed of one Cu-N bond and at least one Cu-O bond. The complete leaching of Cu is further evaluated by the slow-growth method at a time of 2.5 ps, as shown in Figure S12. An energy barrier of 0.84 eV is obtained, indicating that this reaction can occur in mild conditions, but is slightly higher by 0.14 eV than the incomplete leaching of Cu (0.70 eV). Thus, at a dynamic electrochemical interface, Cu should exist in two transient states, i.e., incomplete leaching combined with one N atom and complete leaching dissolved in water. As the applied potential becomes more negative, H adsorption is easier, leading to a more indirect promotion of Cu leaching (see Figure S13). The projected crystal orbital Hamilton population (pCOHP) 42,43 is used to analyze the bond strength between the central copper atom and the nitrogen atom as H adsorption increases. As shown in Figure 3e, the contribution of the antibonding orbital population below the Fermi level increases with the increase in the number of H adsorbed, which weakens the strength of the Cu-N bond. In addition, the integrated pCOHP (IpCOHP) can quantitatively provide bond strength information. 44 The IpCOHP values of one Cu-N bond in the Cu-N 4 -C are given in Figure 3e, corresponding to -1.58 for the bare surface, -0.95 for 1H, and -0.93 eV for 2H. Obviously, as the number of H adsorbed increases, the absolute value of IpCOHP becomes smaller, indicating that the Cu-N bond becomes weaker, which leads to the transition from difficult to spontaneous leaching of the Cu single atom. Therefore, the adsorption of H is a vital driving force that leads to the leaching of the Cu SA from the catalyst surface. The dynamic evolution of the bond length of Cu-N 4(x) (x = 1, 2, 3, and 4) and pCOHP analysis under the conditions of 3H and 4H coadsorption is shown in Figures S14 andS15, respectively. The effect on the Cu leaching behavior is consistent with the results of the two H adsorption. In addition to the adsorption of H on N sites, the adsorbates also have an important effect on the Cu single atom itself. Five intermediates (H, CO, CO 2 , H 2 O, and OH) are considered for adsorption on the Cu single atom (Figure S16). The strength of the weakening of the Cu-N bond is CO > H(Cu) > H(N) > OH > CO 2 > H 2 O. However, due to the lack of saturation of N, the adsorption of CO and H on Cu single atoms (0.79 and 2.03 eV, respectively) is unfavorable compared to H on N (0.28 eV) in the leaching of Cu. Therefore, the adsorption of H on the N site is more likely to promote the leaching of the Cu single atom from the catalytic surface.

As mentioned above, some Cu atoms remain connected to the N atoms and are tethered to the catalyst surface, while some Cu atoms dissolve in the aqueous solution. Under a dynamic environment, the collision of the Cu atoms in the two states forms a transient Cu 3/4 cluster structure as the true catalytic active center. We further simulate the aggregation process of Cu atoms under working conditions. It is observed in Figure 3f that the two Cu atoms aggregate is less than 2.5 ps during AIMD and the Cu (1) atom remains bound to N (1) and N (2) . At 15 ps, the bond length of Cu (1) -Cu (2) was further shortened to 2.26 Å, and the bond lengths of Cu (1) -N (3) and Cu (1) -N (4) were stretched to 3.09 and 3.03 Å, respectively. Due to such a high computational cost, the dynamic simulation calculations on the formation of Cu 3/4 small clusters had to be stopped, but we believe that the Cu 3/4 cluster will be observed in long simulations in a dynamic solid-liquid interface environment. Once more and more Cu small clusters are formed, they can accelerate the CO 2 reduction to ethanol. Our previous study has proved that the Cu 4 cluster supported on g-C 3 N 4 can promote the reduction of CO 2 to ethanol. 45 To date, the restructuring dynamic process of Cu single atoms forming Cu clusters has been clearly understood during the CO 2 reduction reaction.

After the CO 2 reduction reaction is over, the electrode potential is switched off to stop the electrocatalytic reaction or to +1.0 V vs RHE. At this time, the Cu small clusters (Cu 0 ) can be restored to atomically dispersed copper (Cu 2+ ). 38,39,46 Compared to the leaching process of Cu SA, the return of the Cu atom is more complex and challenging in the simulation. The formation of Cu-N 4 sites is accompanied by oxidation of the metallic copper under the positive electrode potential, resulting in the oxidation state of Cu increasing to +2. Here, simple Cu 4 and Cu SA-Cu 3 models are used to describe the electrochemical deposition process in thermodynamics (see Figure S17). The free energy of the Cu SA-Cu 3 formation becomes more negative as the potential increases, which means that the production of Cu SA-Cu 3 is more favorable, as shown in Figure 4a. Recently, an experimental study has confirmed that the occurrence of fast oxygen exchange between HCO 3 and H 2 O facilitates the formation of the highly oxidative hydroxyl radicals (OH • ), which leads to the promotion of the fast reoxidation of Cu. 47 Therefore, we consider the introduction of two OH • radicals in the system to assess the reoxidation process of the Cu 4 cluster at +1.0 V vs RHE. The two configurations, including a Cu 4 cluster combined with two OH • radicals (Cu 4 (OH • ) 2 ) and a Cu 3 cluster with an escaped Cu atom bonded to two OH • radicals (Cu 3 + Cu(OH • ) 2 ), are constructed and performed for AIMD simulations of 6 and 2.5 ps, respectively. Figures 4b andS18 show that two OH • radicals oxidize the first coordination shell of Cu-Cu after 2 and 0.5 ps, respectively. After this, the Cu (1) atom bonded to the catalyst substrate is farther and farther away from the other three or two Cu atoms. For Cu 3 + Cu(OH • ) 2 , the OH • radical runs completely onto the Cu 3 cluster to carry out further oxidation after 0.5 ps, while the escaped Cu is coordinated with the three H 2 O molecules (Cu 3 (OH • ) 2 + Cu(H 2 O) 3 ). Therefore, the presence of OH • radicals plays a dominant role in the fast reoxidation of Cu compared to pure aqueous solutions (Figure S19). It is noted that the net charge of the system is in relative equilibrium during the oxidation process, and almost no additional electrons are transferred to the electrodes, indicating that the Cu cluster is oxidized to Cu δ+ by OH • radicals. After the formation of Cu(H 2 O) 3 , the Cu redeposition process occurs when Cu(H 2 O) 3 diffuses to the adjacency of the N 4 -C site, i.e., Cu migrates from an aqueous solution to be coordinated with a N atom as step #1, overcoming a kinetic energy barrier of 0.85 eV (Figure 4c). The Cu-O (1,2,3) and Cu-N bonds (Figure S20) are 2.12, 2.15, 2.20, and 2.06 Å, respectively. Furthermore, similar to the Cu leaching process, to avoid the complexity of multiple reaction paths, the He atom is again introduced directly below the Cu atom (the Z-direction of the He atom is fixed). The reaction path of the Cu back to the N 4 -C vacancy process is simplified to a shortening of the Cu-He bond as step #2, overcoming an energy barrier of 0.50 eV. During this reaction, three H 2 O molecules return to the aqueous solution with an average bond length of 1.94 Å for Cu-N bonds (Figure S21), consistent with the experimental characterization for the initial Cu-N bond length (1.95 Å). 18 That is, once the Cu atom is coordinated to one N, the Cu atom quickly returns to its initially dispersed state. Therefore, the synergistic effect between highly oxidative OH • and the strong chelating ability of the N 4 -C site under a positive potential promotes the return of the Cu cluster to the Cu SA state, completing the cycle of restructuring behavior. Figure 4d shows a dynamic and reversible transformation mechanism between the Cu SA to the Cu 3/4 small clusters during electrolysis. It is mainly divided into two processes, that is, the escape and return of the Cu SA. The adsorption of H under a negative potential promotes the conversion of the Cu SA to Cu clusters, thereby optimizing the performance of the catalyst and enhancing the selectivity of ethanol. When the electrode potential is switched off or to a positive value, Cu SA returns to the atomic dispersion through oxidation and redeposition processes. To better understand the specificity of the structural transformation of Cu-N 4 -C during CO 2 reduction, we also consider the free energy of one H adsorption on a Cu single atom coordinated with O 4 -C and Fe, Co, and Ni coordinated with N 4 -C at U = -1.2 V. The calculated results indicate that the Cu-N bond is the weakest among all of the M-N bonds, which is consistent with the metal leaching phenomenon reported in experiments for Cu only during CO 2 RR. [18][19][20][21] More details can be found in the Supporting Information (Figures S22 andS23).

CONCLUSIONS

In summary, we have systematically investigated the dynamic stability of the Cu single-atom catalyst under working conditions, namely, the reversible transformation between Cu SA and Cu 3/4 clusters, using the "constant-potential hybridsolvation dynamic model". We reveal that the adsorption of H, depending on the potential, is a crucial driving force for the conversion of Cu single atoms into Cu clusters, while the reversion of Cu clusters to the atomic dispersion state is dominated by oxidation reactions, and OH • plays a crucial role in this process. These results are impossible to be identified using the previous charge-neutral DFT calculations. As a result, considering the effect of surface charge and explicit water solvent is an important means to understand the dynamic stability of SACs under working conditions. We expect a systematic exploration of the stability of SACs in real electrochemical reactions in the future. In addition to SACs, the studies on the oxidation state and surface reconstruction of copper catalysts under CO 2 reduction conditions have been recently reported experimentally. [48][49][50] Although these studies reveal that the applied electrode potential is the dominant factor in the surface reconstruction of copper catalysts during CO 2 reduction, the mechanism for the reconstruction of the catalysts is unclear. Therefore, an in-depth understanding of the dynamic stability of catalysts, whether SACs or bulk materials, can provide not only a fundamental understanding of the evolution of structural and electrical properties under operating conditions but also new ideas for the design and development of catalysts with better performance.

■ ASSOCIATED CONTENT

* sı Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/jacs.2c07178.

Construction of models and computation details; raw data from AIMD and the constrained MD simulations; dynamic evolution of the bond length and surface charge; free energies and activation energies for H adsorption; and pCOHP comparison for xH adsorption at N sites (where x = 0, 1, 2, 3, and 4) (PDF)

■