Sulfur-enhanced dynamics of coinage metal(111) surfaces: Step edges versus terraces as locations for metal-sulfur complex formation

Coarsening or ripening of nonequilibrium nanostructured systems is a generic phenomenon related to the tendency of these systems to reduce the energy penalty associated with a high density of interfaces by increasing the characteristic length scale of structures.¹ One broad class of such systems involves arrays of metal nanoclusters (NCs) either in solution² or supported on surfaces,^3–5 where the above-mentioned energy penalty can be regarded as associated with broken bonds at the periphery of the NCs. In this case, coarsening involves a decrease in the time of the number (or density) of such NCs and a corresponding increase in their mean size, consistent with mass conservation. A classic mechanism of coarsening is Ostwald ripening (OR) wherein mass is transferred from smaller than average to larger than average NCs, which results in the ultimate disappearance of the former.³ A competing mechanism is NC or particle migration and coalescence (PMC), also called Smoluchowski ripening (SR),⁵ although the current contribution is focused on OR type processes.

Further restricting attention to supported metal nanoclusters on flat single-crystal surfaces, the default expectation is that OR under pristine ultrahigh-vacuum conditions will involve the detachment of individual metal atoms from smaller NCs, followed by diffusion across the surface and attachment to larger NCs.³ There is much interest in the understanding of such OR (or SR) processes as they provide mechanisms for the degradation of the catalytic properties of 3D metal NCs supported on oxides or other substrates,⁶ noting that OR reduces the NC surface area and, thus, catalytic activity. However, more fundamental and definitive studies of OR (and SR) have considered 2D metal NCs in metal homoepitaxial systems, where a combination of scanning tunneling microscopy (STM) and high-level atomistic-level modeling can provide a detailed characterization and understanding of coarsening. The reader is referred to related reviews by the Thiel⁵ and Morgenstern⁴ groups.

Applications to catalysis motivate consideration of coarsening processes under nonpristine “operating conditions” where nonmetallic chemisorbed species are present. In this context, we note the existence of a long-standing hypothesis in catalysis community that the degradation of supported 3D metal NC catalysts can be enhanced by the formation of volatile metal-oxygen complexes that facilitate metal mass transport and, thus, enhance coarsening.^7,8 However, rather than these more complex catalytic systems, simpler metal homoepitaxial systems incorporating 2D metal NCs which will be considered here again offer the potential for a more definitive analysis of such complex-facilitated mass transport and coarsening.⁹ In the context of these simpler systems, we also note recent general interest in refining the traditional view that metal surfaces are static or frozen during chemisorption or catalytic surface reaction to instead account for their potentially fluxional nature.¹⁰ The above considerations of how chemisorbed species can enhance surface mass transport and coarsening should be regarded as contributing to this more general issue.

Henceforth, in this contribution, we will focus on the dynamics of metal homoepitaxial systems, which involve a low-index metal surface generally supporting 2D metal NCs and also including extended step edges, in the presence of chemisorbed chalcogens. Both serendipitous and controlled studies have revealed chalcogen-enhanced surface dynamics on coinage metal (M) surfaces, particularly on M(111) surfaces, but to some extent also on M(100) surfaces.⁹ Early observations involved accidental exposure of Au(111) to air by Peale and Cooper¹¹ and of Ag(100) by Layson and Thiel,^12,13 were most likely oxygen-induced enhancement. More specifically, these studies, and those described below, generally explored the rate of decay in size of 2D homoepitaxial metal NCs nearby extended step edges, observing much higher decay rates in the presence of chemisorbed chalcogen relative to the pristine chalcogen-free system. Subsequently, a detailed and controlled study was performed at Sandia National Laboratory by Ling et al.¹⁴ for S on Cu(111) revealing dramatic enhancement of NC decay rates even with trace amounts of S. Supporting density functional theory (DFT) analysis suggested that mass transport was promoted by the formation of Cu₃S₃ “trimer” complexes.¹⁵ This proposal was later questioned with linear CuS₂ and Cu₂S₃ “hearts” being suggested as possible alternative mass transport agents.¹⁶ Similar dramatic enhancement for trace S on Ag(111) was observed by the Thiel group at Iowa State University, where AgS₂, Ag₃S₃, and other complexes were identified as possible as promotors of mass transport.¹⁷ In contrast, the Thiel group observed only moderate enhancement even for larger amounts of S on Au(111), where AuS₂ and Au₄S₄ were suggested as possible mass transport agents.¹⁸ 

The rationalization for enhanced surface mass transport due to metal-sulfur complexes is that these have low formation energy (i.e., high volatility) relative to that of metal adatoms (which must be extracted from steps on the surface).^13,15 Note that formation energies determine equilibrium adspecies concentrations on the surface. Especially, lower complex formation energies lead to higher surface populations, which can more than compensate for the anticipated lower mobility of the complex again relative to that of adatoms, thereby potentially enhancing mass transport.

In fact over the last decade, motivated by observations of enhanced coarsening, extensive studies have been performed by the Thiel group examining various aspects of behavior for S on the $⟨111⟩$ surfaces of coinage metals (M)^19–21 exploring all of M = Cu,^16,22,23 M = Ag,^17,24–28 and M = Au.^18,29,30 One observation from these studies is the particular stability of linear S-M-S motifs, either as isolated MS₂ complexes or as components of larger complexes. We comment on this feature further in this contribution. In addition, the Thiel group studied and compared behavior for $⟨111⟩$ surfaces with that for S on the $⟨100⟩$ surfaces of M = Cu,³¹ M = Ag,^32,33 and M = Au.³⁴ One conclusion from their experimental studies for these systems, aided by extensive supporting DFT analysis, is as follows:³⁵ complex formation is facile and dramatically enhances mass transport on Cu(111)^14,16 and Ag(111),¹⁷ but only moderately on Au(111).¹⁸ In contrast, complex formation is less facile on the $⟨100⟩$ surfaces of Cu and Ag where little enhancement is observed. Complex formation is predicted to be more facile on Au(100),³⁵ but in this case, studies to assess the existence of S-enhanced surface dynamics have not yet been performed. Figure 1 displays pictures of the Thiel group graduate students involved in this series of work, together with that of their advisor.

We also note one unique feature of the Thiel group efforts, absent in studies by other groups, was the undertaking of ultra-low STM studies aiming to directly image by “freezing in” at 5 K, and at low S coverage, the isolated complexes implicated in enhanced mass transport. This effort was performed in collaboration with the groups of Maki Kawai and Yousoo Kim at RIKEN in Japan. The interpretation was facilitated by a detailed comparison of STM images with DFT-generated images. Interestingly, Cu₂S₃ “hearts” and larger Cu_mS_m+1 concatenated chains (but not Cu₃S₃ “trimers”) were observed on Cu(111).¹⁶ Only relatively large complexes were observed on Ag(111), including Ag₁₆S₁₃, although fragments of these (such as linear AgS₂) could be predominant and facilitate mass transport at higher temperature (T).²⁷ Note that Ag₃S₃ “dots” as part of a “row-dot” structure at higher S coverage were also observed. Only S adatoms were observed on Au(111).^29,30

In this contribution, we assess the kinetics of complex formation focusing on the case of S on Cu(111) and considering both formation on terraces and directly at step edges. Relevant background information is that there is an energetic preference for S to be adsorbed at step edges rather than on terraces. With this feature in mind, it is appropriate to distinguish different regimes of the coverage, θ_S, of S: (r1) very low θ_S where neither A nor B type steps are saturated, and just A steps are populated; (r2) low θ_S where A steps are saturated, but not B steps that are partially populated; and (r3) higher θ_S where all steps are saturated, and excess S resides on terraces. Note that an assessment of the degree to which steps are saturated depends not just on θ_S, but also on the quality of the sample/substrate, i.e., on the step density. Only a very low θ_S is needed to saturate steps on a high-quality substrate. Evaluation of the complex formation energy mentioned above depends upon the regime. The first calculation¹⁵ of the formation for Cu₃S₃ applied only for the regime (r3), although this was not explicitly noted in that study.

In Sec. II, we extract chemical potentials for a variety of structures involving adsorbed S that can form at step edges for increasing amounts of S. These, together with the chemical potentials for complexes, lead to corresponding values of the formation energies for various complexes. Our analysis refines previous assessments. Next, in Sec. III, we address the greater challenge of assessing the kinetics of the complex formation process and its relationship to enhanced coarsening. Note that while complex formation can appear a viable candidate for enhanced surface mass transport based on the assessment of surface thermodynamics (i.e., based upon the complex formation energy) and complex mobility, it may not be kinetically viable. First, we explore the kinetics for possible mechanisms of complex formation on terraces, extending previous studies to identify a new viable pathway for CuS₂ formation. Second, we analyze the kinetics of complex formation by detachment directly from S-decorated step edges or structures at step edges, a previously ignored scenario that is shown to be a viable alternative to the complex formation on terraces. This analysis is enabled by nudged-elastic-band determination of the minimum energy path (MEP) for detachment utilizing machine-learning (ML) derived potentials for the metal-sulfur system retaining near ab initio DFT accuracy. Conclusions are provided in Sec. IV.

It is instructive to first introduce the sulfur chemical potential, μ_S, which provides an assessment of the relative thermodynamic preference of various structures which might be formed on the surface by incorporating adsorbed S. This quantity is determined from DFT calculations using the VASP code with the PAW method and the PBE functional, as described in detail in previous work on these S + metal systems.^20,35 We will utilize results from those studies to perform additional analysis here. The periodic computational unit cell includes a Cu substrate (which could be a perfect flat $⟨111⟩$ surface or a stepped surface), together with an adsorbate structure involving n ≥ 1 S atoms and possibly m ≥ 0 additional Cu atoms on the substrate.

For the case m = 0, which applies for an isolated S adatom on the Cu(111) terrace or S decorated step edges, we introduce the quantity

where E(⋯) denotes the total energy of the indicated system and E₀ is the reference energy chosen here as 1/2 E(S_2,gas). For configurations that involve S atoms in a single configuration or local environment, μ_S* corresponds to the S chemical potential, μ_S = μ_S(S). However, for configurations with S in different environments, μ_S* averages over the energetics of these. In this case, an additional step is required to determine the S chemical potential, μ_S. This process is described below and in  Appendix A noting that μ_S corresponds to the differential increase in (free) energy of the system upon increasing the amount of S. For m > 0, one is typically considering adsorbed Cu-S complexes on terraces, and the S chemical potential associated with the complex is given by^16,20,22

with E(…) and E₀ as above and where μ_Cu denotes the energy of a Cu atom in the bulk (i.e., the negative of the bulk cohesive energy plus a Cu atom self-energy).

The highest value of μ_S in the low S coverage regime is roughly μ_S(T) = −1.91 eV corresponding to isolated S adatoms on terraces (T).¹⁶ The lowest value of μ_S is μ_S(A) = μ_S*(A) = −2.24 eV, which is achieved for S at $⟨100⟩$ microfaceted A-type steps²² up to the point where such steps are saturated by a double-spaced row of S as shown in Fig. 2(a). Thus, μ_S = −2.24 eV for regime r1 with unsaturated A steps. For higher S coverages where A-type steps are fully saturated, and $⟨111⟩$ microfaceted B-type steps begin to be populated again by a double-spaced row up to saturation, i.e., regime r2, μ_S is given by μ_S(B) ≈ −2.12 eV.²² In this regime r2, one finds essentially identical μ_S for a double-spaced row of S forming either along the step edge or just above the step edge, as shown in Figs. 2(b) and 2(c). States in regime r2 with a mixture of double-spaced S at and just above the step should also have essentially the same μ_S ≈ −2.12 eV.

For somewhat higher amounts of S at the onset of regime r3, pairs of rows of S at and just above step edges form, sometimes inducing local reconstruction at the step edge.²² The latter allows S to reside at preferred fourfold hollow (4fh) sites³⁶ (see Fig. 3). Here, we determine not just μ_S* (which averages over S in different environments) but also the chemical potential, μ_S, which should be regarded as associated with the additional S atoms beyond those shown in Fig. 2. For example for the A step, one has that μ_S* = 1/2 [μ_S(A) + μ_S], where μ_S is associated with the additional S atoms in the row just above the step. See  Appendix A for an alternative derivation. Apart from the case in Fig. 3(a), these μ_S are still significantly lower than μ_S(T), so these configurations are preferred versus populating S adatoms on terraces.

Next, consider adding further S sufficient to form triples of rows with one at and two just above the step edge (not shown). The μ_S associated with the additional S in the third row is not significantly lower than μ_S(T), so dispersion of the additional S on terraces is at least as favorable as these structures. However, an alternative configuration for this amount of S on the surface is to form an extended Cu_mS_m+1 chain just above an A step²² [see Fig. 4(a)]. In this case, we can use the relation μ_S* = 1/3 [μ_S(A) + 2 μ_S] to obtain the chemical potential, μ_S, associated with S in the double row decorating the Cu chain. In addition, complex reconstructions with extended triangular motifs occur at B steps generating multiple 4fh sites preferred by S²² [see Fig. 4(b)]. We will not discuss the latter configurations further here.

The above analysis involves comparing similar chemical potentials, where it should be noted that there is some uncertainty in DFT values. However, a general conclusion is that as more S is added to the surface (but remaining in a low-coverage regime where the S population on terraces is low), both μ_S* and μ_S increase toward μ_S(T) = −1.91 eV. It is also appropriate to note that STM studies by the Thiel group²² have directly imaged a variety of the structures described above, including the S decorated A step in Fig. 2(a) as well as the more exotic chains and reconstructions in Fig. 4.

Finally, in Fig. 5, we present the S chemical potentials, μ_S(Cu_mS_n), for key complexes, obtained from the analysis in Ref. 20 and also show the configurations of these complexes.

First, we consider the scenario where the Cu_mS_n complex is formed by extracting Cu from the “bulk” (in practice from kink sites at step edges) and where the source of S will depend on the coverage regime. The formation energy for the complex on the terrace, E_form(Cu_mS_n), corresponds to the energy difference between the final state with the complex on the terrace and the initial state with Cu in the bulk and S in the relevant source configuration. If δE(Cu_mS_n) = E(Cu_mS_n + substrate) − E(substrate) denotes the total energy associated with the complex, then this formation energy satisfies¹⁶ 

where again μ_Cu denotes the energy of a Cu atom in the bulk, and we note that μ_S = δE(S) − E₀. We emphasize that the appropriate μ_S will depend on the source of the S and, thus, on the S coverage. We note that another formulation of E_form(Cu_mS_n) makes explicit that this quantity is determined by the balance between the energy increase due to forming isolated Cu and S atoms on the terrace and the energy decrease due to binding within the complex (see  Appendix B). Previous analyses have always assigned μ_S = μ_S(T) = −1.91 eV (or a refined value for finite S coverage on terraces), but this formulation only applies in regime r3 with saturated steps and excess S on terraces.

For regime r1 with unsaturated A steps, one should select μ_S = μ_S(A) = −2.24 eV. For regime r2 with saturated A steps but unsaturated B steps, one selects μ_S = μ_S(B) ≈ −2.12 eV. However, in regime r3 with saturated steps and excess S on terraces, one selects μ_S = μ_S(T) ≈ −1.91 eV. Values of E_form(Cu_mS_n) for these three regimes are summarized in Table I. For regime r3, these E_form(Cu_mS_n) differ somewhat from those reported in Ref. 20 where μ_S was selected to incorporate an S coverage dependence versus the assignment here of μ_S = μ_S(T) = −1.91 eV. Table I also reports values for diffusion barriers, E_d, of complexes which will be relevant for the analysis of complex-mediated coarsening in Sec. III.

Although likely not a significant component of enhanced coarsening, one can also assess the formation of small complexes such as CuS₂ or Cu₂S₃ via detachment or extraction from extended Cu_mS_m+1 chains (with large m) above A steps which are shown in Fig. 4(a). Recall that these chains were observed in STM experiments. Contrasting the above analysis, now the source of Cu is not the bulk substrate. Results from such an analysis are presented in  Appendix B.

As indicated in Sec. I, Ostwald ripening (OR) of a distribution of 2D monolayer Cu NCs or islands on Cu(111) involves the transfer of Cu from smaller than average to larger than average islands. Thus, OR reduces the number of islands and increases their mean size, thus reducing the total island perimeter length, and thereby reducing the energy cost associated with “broken bonds” at their perimeters. The coarsening rate associated with mass carriers, C, reflects the product of the mobility of C described by a diffusion coefficient, D_C ∼ exp[−E_d(C)/(k_BT)], where E_d(C) denotes the terrace diffusion barrier for the carrier, and the equilibrium concentration or coverage of C satisfying θ_eq(C) ∼ exp[−E_form(C)/(k_BT)], where E_form(C) denotes the carrier formation energy. Consequently, the effective barrier for OR has the form^17,18,37

We include a possible additional contribution, Δ, to account for the scenario where there is a substantial additional barrier (beyond the terrace diffusion barrier) for attachment of carriers to island edges. However, this term boosting E_OR only appears if the so-called attachment length L_attach = exp[Δ/(k_BT)] (measured in lattice constants) exceeds the mean island separation, L_isl. We also note that the barrier for detachment of a mass carrier, C, from step edges, E_detach, corresponds to E_OR. For a pristine Cu(111) surface, i.e., in the absence of S, the carriers are Cu adatoms, C = Cu, and one has that E_d(Cu) = 0.05 eV and E_form(Cu) = 0.78 eV,¹⁶ and there is no attachment barrier (Δ = 0), so that E_OR(Cu) = 0.83 eV.

For a Cu_mS_n complex to enhance OR, a necessary condition is that E_OR(C) be below E_OR(Cu) = 0.83 eV. Based on this requirement, we can immediately rule out CuS and Cu₂S₂ in regime r3 based on the results in Table I. In regime r2, only CuS₂ can potentially enhance OR. In regime r1, no complexes can enhance OR. However, for complex formation on terraces (where Δ = 0), E_OR(C) < E_OR(Cu) is not sufficient for the complex to enhance OR, and the analysis of kinetics as presented below is required. For complex formation at step edges, E_OR(C) < E_OR(Cu) should be sufficient to enhance OR. However, the challenge is that now generally Δ > 0, and a nontrivial analysis of the kinetics of attachment/detachment is required.

As noted above, the picture presented in the Sandia analysis of complex formation by extraction of Cu atoms from steps and reaction with S on terraces¹⁴ is only viable in regime r3 with saturated steps and excess S on terraces (or in regime r2 for CuS₂). Kinetics is assessed by the analysis of appropriate reaction-diffusion equations (RDEs) for Cu adatom and complex concentrations, where these equations include terms that describe the complex formation, terrace diffusion, and dissolution processes.^5,14,17 The boundary conditions for these equations specify that the concentration of Cu adatoms matches its equilibrium value at NC edges (and thus is higher at smaller NC) and also impose a zero-flux boundary condition for complex concentrations at island edges (since complexes are not formed or removed at step edges).^14,17 The original Sandia study¹⁴ just utilized a simplified coupled linear pair of RDEs for the concentration of Cu atoms and a single complex, C. However, this particularly insightful and elegant analysis did effectively highlight the importance of a “reaction length” for L_R ∼ (D_Cu/k)^1/2, where D_Cu is the Cu adatom diffusion coefficient (which is far higher than that for complexes) and k is a first-order reaction rate for complex formation. If L_R > L_isl (the mean island separation), then the formation of complexes is sufficiently inhibited that the quasi-uniform equilibrium concentration of complexes does not sufficiently strongly couple to gradients in the Cu adatom density to impact mass transport. It is necessary that L_R < L_isl for enhanced transport and coarsening.

Going beyond the above simplified generic analysis, a detailed and predictive system-specific treatment of the kinetics must be based on the appropriate coupled nonlinear RDEs with boundary conditions as described above. Linearization of these RDEs about spatially uniform equilibrium concentrations for various species leads to the natural emergence of reaction lengths together with explicit expressions for these lengths.^5,17,37 A previous analysis following the STM observation of Cu₂S₃ complexes explored the possibility that C = Cu₂S₃, formed by the reaction Cu + CuS₃ → Cu₂S₃, could enhance mass transport.¹⁶ However, the corresponding reaction length far exceeds L_isl, so this does not appear to be a viable mechanism for enhanced OR.

Thus, we consider an alternative pathway for enhanced OR in which CuS₂ dominates mass transport and where CuS₂ is formed on terraces via the trimolecular reaction, Cu + S + S → CuS₂. In this case, the appropriate nonlinear RDE for the surface concentrations or coverages (θ) of Cu adatoms and CuS₂ are

and

Here, F(Cu + S + S) denotes the rate for diffusion-mediated formation of CuS₂ and R(Cu + S + S) denotes the rate of dissociation of CuS₂. These are given by

and

where E_bind(CuS₂) > 0 is the binding energy strength for adsorbed CuS₂ (i.e., the difference in energy between Cu and S dispersed on the terrace and incorporated in CuS₂). This rate selection appropriately satisfies detailed balance. Given that in regime r3, one has

(see  Appendix B), equating F() and R() correctly recovers the equilibrium concentration $θCuS2$(equil) = exp[−E_form(CuS₂)/(k_BT)] (θ_S)². Other possible reaction paths leading to the formation and removal of Cu and CuS₂ are left implicit in (5) and (6).

As indicated above, these nonlinear equations are linearized for small deviations about spatially uniform equilibrium concentrations. Setting θ_Cu = θ_Cu(equil) + δθ_Cu and $θCuS2$ = $θCuS2$(equil) + δ$θCuS2$⁠, from (5) one obtains the linearized equation

with k(Cu + S + S) = (D_Cu + D_S) (θ_S)². As noted above, the activation barrier for terrace diffusion of Cu is E_d(Cu) = 0.05 eV. Our separate calculations for S adatoms show that E_d(S) = 0.15 eV (from DFT using the PBE functional with a 2 × 2 lateral supercell and the Cu{111} surface represented by a three-layer slab), in reasonable agreement with a previous estimate³⁸ of E_d(S) = 0.17 eV. Thus, one has that D_Cu >> D_S and consequently that the relevant reaction length is given by $LR=(DCu/k)1/2≈$ 1/θ_S (in lattice constants for θ_S in monolayers). For θ_S = (5–20) × 10⁻³ monolayers, L_R is well below L_isl in experiments.⁵ The feature that (10) has solutions for δθ_Cu vary exponentially in space with decay length L_R < L_isl demonstrates that there is a sufficiently strong coupling of the Cu adatom concentration to that for CuS₂ complexes to induce complex-enhanced Cu mass transport. Thus, enhancement via this mechanism is operative for both regime r3 where E_OR(CuS₂) = 0.19 eV and for regime r2 where E_OR(CuS₂) = 0.61 eV (vs E_OR(Cu) = 0.83 eV).

There has been no previous analysis of the possible direct formation of complexes at step edges, which specifically aims to assess whether such a mechanism can lead to enhanced OR. The challenge is to identify viable pathways for formation and detachment in what is a vast multi-dimensional phase space of possible configurations at step edges and to determine the MEP and, thus, the barrier for such detachment processes. An additional complication is the large number of Cu–sulfur structures and reconstructions that can occur at step edges and that could provide potential sources for complexes at least in regime r3.

First, we observe that such direct formation of complexes, C, as a pathway for enhanced OR is not possible in regime r1 since E_OR(C) is no smaller than E_d(C) + E_form(C) for which values exceed E_OR(Cu) for all complexes C. Such enhancement is plausible in regime r2 for C = CuS₂ but perhaps unlikely given the possibility of a significant extra barrier for detachment, Δ. Thus, we focus on regime r3 where the higher S chemical potential and the above-mentioned diversity of Cu–sulfur step edge structures, as well as the significant population of S on terraces, should facilitate the formation of various complexes, C, with E_OR(C) < E_OR(Cu).

Given the vast phase space of possible complex formation and detachment pathways, comprehensive exploration at the level of DFT determining MEPs with low barriers is inefficient due to the computational demand. However, there have been significant recent advances in the use of ML strategies to develop force fields retaining high-level near-DFT accuracy.^39,40 In implementing this approach, we have performed ab initio molecular dynamics (AIMD) simulations of the dynamics of Cu nanoparticles with 60–1000 Cu atoms and with various amounts of adsorbed S and also of copper sulfide nanoparticles. The S + Cu(111) system was also simulated mostly for validation and it only contributed to a small part of the training set data. These data, which contain information about how energies and forces depend on atomic coordinates, are then used to train a set of ML potentials. We utilize the DeepMD-tool kit⁴¹ due to its ease-of-use and efficiency. With these potentials, e.g., MD simulation can be performed for much longer times than in AIMD. In our case, much more efficient and extensive nudge elastic band (NEB) analysis⁴² can be performed to determine MEPs for complex formation and detachment. In a separate contribution, we will describe in more detail the development of these potentials and also present extensive MD simulations to probe dynamics at step edges in the S + Cu(111) system.

Here, we just perform a limited but targeted analysis to provide a proof-of-principle demonstration that complexes can be formed directly at or near step edges with a sufficiently low barrier to facilitate enhanced OR. We focus on two cases: (1) detachment of CuS₂ from kink sites at an S-decorated B step and (2) detachment of CuS₂ from extended Cu_mS_m+1 chains just above A steps.

For formation and detachment of complexes at step edges, one naturally anticipates that this process is facilitated at kink sites where the Cu has lower coordination and, thus, is more readily extracted. Thus, in Fig. 6, we show one such example of CuS₂ formation at a kink site on an S-decorated B step, which is facilitated by an additional pair of nearby S adatoms on the terrace. As we are considering regime r3, the population of such terrace S adatoms is nonnegligible. The detachment pathway involves the relatively mobile Cu atom starting at the kink site being passed between relative immobile nearby S adatoms to eventually form a CuS₂ complex on the terrace by inserting itself between the two terrace S adatoms. For this process, the detachment barrier is only E_detach(CuS₂) = 0.43 eV as shown in Fig. 7, well below E_OR(Cu), so it should facilitate enhanced OR. It is interesting to note the formation of a CuS₃-type motif around the transition state of the detachment process. Note that as is appropriate, the difference between the final and initial energies of the MEP roughly corresponds to the DFT estimate of the formation energy for CuS₂ in regime r3. A couple of additional comments are appropriate. First, the energetics of detachment is primarily controlled by the two S atoms at the kink site and the two terrace S atoms and, thus, does not change dramatically if the other 3 S adatoms in Fig. 7 are absent. Second, our analysis of the analogous detachment process at a kink site on an A step yields a significantly higher barrier. This is perhaps not surprising as S adatoms are more stable at the A step than at the B step.

Next, we consider the formation and detachment of complexes by extraction from extended $CumSm+1$ chain above A steps aided by terrace S adatoms, where we explore only the case of CuS₂ formation. One might naturally first consider detachment from the end of the chain. However, an effective NEB analysis requires judicious initial appropriate placement of the terrace S adatoms, and an appropriate selection (of the many possible choices) is not obvious. Thus, here we just consider the process of extraction CuS₂ from the middle of an extended $CumSm+1$ chain where an initial appropriate placement of the terrace S adatoms is more clear. Figure 8 shows one such example of this process facilitated by a pair of nearby S adatoms on the terrace, as is reasonable in regime r3. Analogous to detachment from a kink site at a step edge, the pathway involves passing a relatively mobile Cu atom between relative immobile nearby S adatoms. The MEP shown in Fig. 9 suggests a detachment barrier of ΔE(CuS₂) = 0.48 eV. However, our separate analysis based on the same ML generated potential reveals that there is significant repulsion of about E_rep = 0.11 eV between the terrace and chain S in the initial configuration. Thus, starting with terrace S further from the chain would lower the initial energy on the MEP by E_rep and boost the detachment barrier to E_detach(CuS₂) = E_rep + ΔE = 0.59 eV, still significantly below E_OR(Cu).

A detailed analysis of coarsening kinetics for the mechanism involving direct detachment of complexes from step edges of nanoclusters would be based upon conventional diffusion equations for complex concentrations [analogous to classic OR treatment that applies for coarsening of islands on S-free metal(111) surfaces]. In contrast to the S-free case where the Cu concentration matches its equilibrium value at step edges, now the boundary conditions on the complex concentration reflect the presence of any additional attachment barriers. Thus, the process can be attachment-limited versus diffusion-limited (the latter applying when Cu adatoms are mass carriers). The boundary conditions should also reflect the feature that the Cu chemical potential is slightly higher at the edge of smaller islands, and this should induce slightly lower detachment barriers than for larger islands. One other complication is that there can be numerous energetically similar detachment pathways involving different arrangements of S adatoms at the step edge and nearby on the terrace. Thus, the determination of the prefactor for the detachment rate, which could exceed typical values, is nontrivial.

For complex formation via detachment from linear Cu_mS_m+1 chains, it is less clear whether this process can effectively contribute to coarsening since the chemical potential of Cu in such chains should be independent of whether they are adjacent to the step edges of large or small Cu islands. If the S concentration is slightly higher at smaller islands, this feature could boost the detachment rate.

Coarsening of supported metal nanoclusters is a phenomenon of fundamental interest,^3–5 but also of importance for applications including the understanding of catalyst degradation.⁶ Of particular interest is the effect of nonmetallic chemisorbed species that are likely present for most applications under “operating conditions.” In particular, it has been suggested that coarsening and more generally surface mass transport can be enhanced due to the formation of adsorbed complexes incorporating these chemisorbed species and metal atoms. This idea has been proposed for complex catalytic systems involving 3D metal nanoclusters on nonmetallic supports⁷ as well as for chemisorption onto low-index metal surfaces supporting 2D homoepitaxial nanoclusters.^9,15 The latter class of systems provides the opportunity for the development of a fundamental understanding and characterization of such enhanced coarsening processes. For S chemisorbed on coinage metal surfaces, following a study by the Sandia group for S/Cu(111),¹⁴ the Thiel group together with theoretical collaborators undertook a systematic study of this entire class of systems. A general goal was to provide a broad set of observations of behavior in different systems, which could underpin a general understanding of requirements for enhanced coarsening. A specific goal was to provide definitive experimental imaging of adsorbed metal-sulfur complexes to validate their existence at least in some systems. These goals were achieved in an extensive series of publications,^16–36 providing a sound framework for the analysis of enhanced coarsening in these systems.

The default paradigm for enhanced coarsening for S on metal surfaces is that complex formation occurs by the extraction of metal atoms from the edges of nanoclusters and their reaction with excess S adatoms on terraces.¹⁴ While this is often plausible from a thermodynamic perspective, an appropriate analysis of the kinetics had not yet demonstrated viability for any specific system. However, such an analysis was achieved here for S on Cu(111), specifically demonstrating that coarsening could be enhanced by the formation of CuS₂ via a trimolecular reaction Cu + S + S → CuS₂. In addition, we introduce a new paradigm of the direct formation of complexes at the step edges of nanoclusters. We demonstrate that this is a viable mechanism for enhanced coarsening. Aided by reliable ML-derived potentials for the Cu-S system, we perform an NEB analysis of the MEP for the formation and detachment process for CuS₂ aided by nearby terrace S. This analysis does find scenarios where the barrier for this process is not prohibitive, at least when there is a significant population of terrace S.

We have noted in the Introduction that our previous studies suggested the particular stability of linear S-Cu-S motifs, either as isolated CuS₂ complexes or as components of larger complexes or chains. The feature that CuS₂ is identified in this study as a natural candidate to enhance OR is compatible with this view. An alternative supporting perspective on the relevant energetics of CuS₂ and other Cu_mS_n complexes is provided in  Appendix C utilizing a many-body cluster expansion for lateral interactions.⁴³ 

Finally, we remark that for the observed enhanced coarsening in S + Ag(111),¹⁷ a similar picture to that for S + Cu(111) is plausible. The mechanism Ag + S + S → AgS₂ on terraces should be effective in enhancing OR. Plausibly, complex formation via direct detachment from S-decorated step edges can also enhance OR.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Division of Chemical, Sciences, Geosciences, and Biological Sciences at Ames Laboratory under the Computational and Theoretical Chemistry (CTC) project. It was performed at Ames Laboratory, which is operated by Iowa State University under Contract No. DE-07CH11358. The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and from the corresponding authors upon reasonable request.

The S chemical potential satisfies μ_S = ∂F/∂N_S, where F is the free energy of the system (approximated here as the total energy, E) and N_S is the number of S adatoms. For A steps decorated by one double-spaced row of S adatoms (regime r1), we neglect interactions between the S. Then, the system energy will vary linearly with the number of S at the step, and for convenience, μ_S is determined from the configuration in Fig. 2(a) where the A step is fully populated (i.e., saturated) by S. Increasing the amount of S to reach regime r2, again the system energy varies linearly with the number of additional S incorporated at B steps to form a double-spaced row. The rate of variation (i.e., the chemical potential) is different than in regime r1 due to the different local environments of S adatoms. For convenience, μ_S is determined from the configuration in Fig. 2(b) or 2(c) where the B step is fully populated by a double-spaced row of S either at or above the step.

At the onset of regime r3 for the scenario shown in Fig. 3(b) or 3(c), the second row of double-spaced S adatoms begins to form just above the A or B steps. Consider determination of μ_S associated with the population of the row of S adatoms just above the reconstructed A step. In Sec. II, we have used the configuration in Fig. 3(b) with a complete double-spaced row both at and above the step to evaluate the chemical potential μ_S from μ_S* = ½ [μ_S(A) + μ_S]. A derivation of this result based upon the relation μ_S = ∂F/∂N_S follows. Consider an A step of fixed length, L_A, in lattice constants, which is fully populated by a double-spaced row of N_S⁻ S adatoms at the step and partially populated by a shorter double-spaced row of N_S⁺ S adatoms just above the step. Then, one has N_S = N_S⁺ + N_S⁻ with N_S⁺ < N_S⁻, where N_S⁻ = ½ L_A is fixed. This configuration has (i) a portion with a pair of double-spaced rows each with N_S⁺ S atoms at and above the step edge associated with mean chemical potential μ_S* and (ii) a portion with a single double-spaced row of N_S⁻−N_S⁺ S atoms at the step edge associated with chemical potential μ_S(A). Thus, it follows that

which recovers the above relation involving μ_S noting that

A similar analysis applies for μ_S in other configurations such as that in Fig. 4(a).

First, consider the formation of complexes on terraces where Cu atoms are extracted from the bulk. We define δE(X) = E(X + substrate)−E(substrate) corresponding to the total energy associated with an adsorbed species, X, which includes the self-energy of the constituent atoms. Then, the total “surface atomization” or binding energy strength E_bind(Cu_mS_n) > 0 of Cu and S atoms within the complex is given by²⁰ 

where Cu_ads and S_ads denote isolated adsorbed atoms on terraces. Formation energies for Cu and S adatoms on terraces (where Cu is extracted from the bulk) are given by

respectively. Here, μ_Cu denotes the total energy of bulk Cu, and δE(S) depends on the source of S. For example, E_form(S) = 0 where the source of S is on the terrace, but E_form(S) > 0 for which the source of S is at unsaturated steps. Then, since the formation energy for the Cu_mS_n complex on terraces by extraction of Cu from the bulk satisfies E_form(Cu_mS_n) = δE(Cu_mS_n) − m μ_Cu− n δE(S), it follows that

Thus, E_form(Cu_mS_n) is determined by the balance between the energy decrease due to binding within the complex and the energy increase due to forming isolated Cu and S atoms on the terrace.

Second, consider the formation of small complexes such as CuS₂ and Cu₂S₃ via extraction from an extended Cu_mS_m+1 chain for large m. From the definition of the S chemical potential, μ_S(Cu_mS_m+1) = μ_S(chain) for the extended chain, which will be independent of large m, it follows from Eq. (2) that

We will also use that δE(S) = μ_S + E₀ and select μ_S = μ_S(T) ≈ −1.91 eV and μ_S(chain) = −1.85 eV below. For extraction of CuS₂ from the end of the chain, via the mechanism Cu_mS_m+1 + S → Cu_m−1S_m + CuS₂, the formation energy satisfies

For extraction from the interior of the chain via the mechanism Cu_m+m′S_m+m′+1 + 2S → Cu_m−1S_m + CuS₂ + Cu_m′S_m′+1, the formation satisfies

A perspective on the energetics of Cu_mS_n complexes follows from the application of a many-body cluster expansion commonly used in lattice-gas models for chemisorbed species.⁴³ Given the presumed particular stability of linear S-Cu-S motifs, it is natural to decompose the total energy of complexes as the sum over (i) an attractive pair interaction of strength ϕ_CuS > 0 for each neighboring Cu–S pair; (ii) an attractive trio interaction of strength ϕ_SCuS > 0 for each linear S-Cu-S unit; and (iii) an attractive contribution of strength ϕ_CuCu > 0 for each nearest-neighbor (NN) pair of Cu adatoms. Within this model, binding energy strengths are given by

To obtain DFT values of these quantities, we use the relation E_bind(Cu_mS_n) = −E_form(Cu_mS_n) + m E_form(Cu), with the DFT value of E_form(Cu) = 0.78 eV. Then, from Table I for regime r3, we obtain E_bind(CuS) = 0.12 eV, E_bind(CuS₂) = 0.64 eV, and E_bind(Cu₂S₃) = 1.53 eV. Then, (C1) provides three equations for the three unknown interaction parameters, the solution of which yields

These values are reasonable, the large value of ϕ_SCuS being consistent with the perceived stability of linear S-Cu-S motifs. Also, the value of ϕ_CuCu reasonably consistent with the DFT value of E_form(Cu), which should correspond to 3ϕ_CuCu.

With these energetic parameters, the model predicts a large value for

reflecting in part the multiple S-Cu-S motifs in Cu₃S₃ (cf. the DFT value of 2.13 eV) and much lower values for

noting the lack of linear S-Cu-S motifs (cf. DFT values of 0.56 and 0.58 eV, respectively).