The propensity of trace amounts of sulfur adsorbed on coinage metal(111) surfaces to dramatically enhance surface dynamics has been demonstrated by STM observations of accelerated 2D island decay for Cu and Ag. It is generally accepted that this enhancement is due to the formation of adsorbed metal-sulfur complexes, which facilitate surface mass transport of the metal. These complexes were originally proposed to form on terraces following the extraction of metal atoms from step edges and subsequent combination with sulfur on the terraces. However, even when thermodynamically feasible, this mechanism may not be kinetically viable for some complexes due to limited coupling of the complex concentration to the surface diffusion flux of metal atoms. Focusing on the case of Cu, we assess various scenarios where complexes are formed either on terraces or instead directly at step edges, the latter being a new paradigm. A new pathway is proposed for the formation on terraces. A rich variety of structures incorporating S at step edges exist, which could provide a viable source for complexes, at least from a thermodynamic perspective. However, it is necessary to also assess the activation barrier for complex formation and detachment from step edges. This is facilitated by the nudged-elastic-band analysis of the minimum energy path for this process utilizing machine-learning derived potentials based on density functional theory energetics for the metal-sulfur system.

## I. INTRODUCTION

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.^{1} One broad class of such systems involves arrays of metal nanoclusters (NCs) either in solution^{2} 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.^{3} A competing mechanism is NC or particle migration and coalescence (PMC), also called Smoluchowski ripening (SR),^{5} 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.^{3} 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,^{6} 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^{5} and Morgenstern^{4} 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.^{9} 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.^{10} 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.^{9} Early observations involved accidental exposure of Au(111) to air by Peale and Cooper^{11} 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.*^{14} 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_{3}S_{3} “trimer” complexes.^{15} This proposal was later questioned with linear CuS_{2} and Cu_{2}S_{3} “hearts” being suggested as possible alternative mass transport agents.^{16} Similar dramatic enhancement for trace S on Ag(111) was observed by the Thiel group at Iowa State University, where AgS_{2}, Ag_{3}S_{3}, and other complexes were identified as possible as promotors of mass transport.^{17} In contrast, the Thiel group observed only moderate enhancement even for larger amounts of S on Au(111), where AuS_{2} and Au_{4}S_{4} were suggested as possible mass transport agents.^{18}

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 $\u27e8111\u27e9$ 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_{2} 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 $\u27e8111\u27e9$ surfaces with that for S on the $\u27e8100\u27e9$ surfaces of M = Cu,^{31} M = Ag,^{32,33} and M = Au.^{34} One conclusion from their experimental studies for these systems, aided by extensive supporting DFT analysis, is as follows:^{35} complex formation is facile and dramatically enhances mass transport on Cu(111)^{14,16} and Ag(111),^{17} but only moderately on Au(111).^{18} In contrast, complex formation is less facile on the $\u27e8100\u27e9$ surfaces of Cu and Ag where little enhancement is observed. Complex formation is predicted to be more facile on Au(100),^{35} 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_{2}S_{3} “hearts” and larger Cu_{m}S_{m+1} concatenated chains (but not Cu_{3}S_{3} “trimers”) were observed on Cu(111).^{16} Only relatively large complexes were observed on Ag(111), including Ag_{16}S_{13}, although fragments of these (such as linear AgS_{2}) could be predominant and facilitate mass transport at higher temperature (T).^{27} Note that Ag_{3}S_{3} “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^{15} of the formation for Cu_{3}S_{3} 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_{2} 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.

## II. ENERGETICS FOR STEP EDGE STRUCTURES AND M_{m}S_{n} COMPLEXES

### A. Sulfur chemical potentials for step edge structures and complexes

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 $\u27e8111\u27e9$ 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_{0} 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_{0} 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).^{16} The *lowest value* of *μ*_{S} is *μ*_{S}(A) = *μ*_{S}*(A) = −2.24 eV, which is achieved for S at $\u27e8100\u27e9$ microfaceted A-type steps^{22} 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 $\u27e8111\u27e9$ 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.^{22} 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.^{22} The latter allows S to reside at preferred fourfold hollow (4fh) sites^{36} (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_{m}S_{m+1} chain just above an A step^{22} [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^{22} [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^{22} 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.

### B. Complex formation energies

First, we consider the scenario where the Cu_{m}S_{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_{m}S_{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_{m}S_{n}) = E(Cu_{m}S_{n} + substrate) − E(substrate) denotes the total energy associated with the complex, then this formation energy satisfies^{16}

where again *μ*_{Cu} denotes the energy of a Cu atom in the bulk, and we note that *μ*_{S} = δE(S) − E_{0}. 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_{m}S_{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_{m}S_{n}) for these three regimes are summarized in Table I. For regime r3, these E_{form}(Cu_{m}S_{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.

. | CuS . | CuS_{2}
. | Cu_{2}S_{2}
. | CuS_{3}
. | Cu_{2}S_{3}
. | Cu_{3}S_{3}
. |
---|---|---|---|---|---|---|

E_{form} regime r1 | +0.99 | +0.80 | +1.64 | +1.20 | +1.02 | +1.20 |

E_{form} regime r2 | +0.87 | +0.56 | +1.40 | +0.84 | +0.66 | +0.84 |

E_{form} regime r3 | +0.66 | +0.14 | +0.98 | +0.22 | +0.03 | +0.21 |

E_{d}(Cu_{m}S_{n}) | +0.33 | +0.05 | — | +0.36 | +0.35 | +0.36 |

. | CuS . | CuS_{2}
. | Cu_{2}S_{2}
. | CuS_{3}
. | Cu_{2}S_{3}
. | Cu_{3}S_{3}
. |
---|---|---|---|---|---|---|

E_{form} regime r1 | +0.99 | +0.80 | +1.64 | +1.20 | +1.02 | +1.20 |

E_{form} regime r2 | +0.87 | +0.56 | +1.40 | +0.84 | +0.66 | +0.84 |

E_{form} regime r3 | +0.66 | +0.14 | +0.98 | +0.22 | +0.03 | +0.21 |

E_{d}(Cu_{m}S_{n}) | +0.33 | +0.05 | — | +0.36 | +0.35 | +0.36 |

Although likely not a significant component of enhanced coarsening, one can also assess the formation of small complexes such as CuS_{2} or Cu_{2}S_{3} via detachment or extraction from extended Cu_{m}S_{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.

## III. ANALYSIS OF COARSENING KINETICS FOR S + Cu/Cu(111)

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_{B}T)], 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_{B}T)], 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_{B}T)] (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,^{16} and there is no attachment barrier (Δ = 0), so that E_{OR}(Cu) = 0.83 eV.

For a Cu_{m}S_{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_{2}S_{2} in regime r3 based on the results in Table I. In regime r2, only CuS_{2} 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.

### A. Complex formation on terraces

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^{14} is only viable in regime r3 with saturated steps and excess S on terraces (or in regime r2 for CuS_{2}). 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^{14} 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_{2}S_{3} complexes explored the possibility that C = Cu_{2}S_{3}, formed by the reaction Cu + CuS_{3} → Cu_{2}S_{3}, could enhance mass transport.^{16} 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_{2} dominates mass transport and where CuS_{2} is formed on terraces via the trimolecular reaction, Cu + S + S → CuS_{2}. In this case, the appropriate nonlinear RDE for the surface concentrations or coverages (*θ*) of Cu adatoms and CuS_{2} are

and

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

and

where E_{bind}(CuS_{2}) > 0 is the binding energy strength for adsorbed CuS_{2} (i.e., the difference in energy between Cu and S dispersed on the terrace and incorporated in CuS_{2}). 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 $\theta CuS2$(equil) = exp[−E_{form}(CuS_{2})/(k_{B}T)] (*θ*_{S})^{2}. Other possible reaction paths leading to the formation and removal of Cu and CuS_{2} 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 $\theta CuS2$ = $\theta CuS2$(equil) + δ$\theta CuS2$, from (5) one obtains the linearized equation

with k(Cu + S + S) = (D_{Cu} + D_{S}) (*θ*_{S})^{2}. 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^{38} 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\u2248$ 1/*θ*_{S} (in lattice constants for *θ*_{S} in monolayers). For *θ*_{S} = (5–20) × 10^{−3} monolayers, L_{R} is well below L_{isl} in experiments.^{5} 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_{2} complexes to induce complex-enhanced Cu mass transport. Thus, enhancement via this mechanism is operative for both regime r3 where E_{OR}(CuS_{2}) = 0.19 eV and for regime r2 where E_{OR}(CuS_{2}) = 0.61 eV (vs E_{OR}(Cu) = 0.83 eV).

### B. Complex formation directly at steps

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_{2} 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^{41} 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^{42} 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_{2} from kink sites at an S-decorated B step and (2) detachment of CuS_{2} from extended Cu_{m}S_{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_{2} 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_{2} complex on the terrace by inserting itself between the two terrace S adatoms. For this process, the detachment barrier is only E_{detach}(CuS_{2}) = 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_{3}-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_{2} 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_{2} 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_{2} 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_{2}) = 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_{2}) = 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_{m}S_{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.

## IV. CONCLUSIONS

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.^{6} 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^{7} 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),^{14} 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.^{14} 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_{2} via a trimolecular reaction Cu + S + S → CuS_{2}. 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_{2} 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_{2} complexes or as components of larger complexes or chains. The feature that CuS_{2} 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_{2} and other Cu_{m}S_{n} complexes is provided in Appendix C utilizing a many-body cluster expansion for lateral interactions.^{43}

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

### APPENDIX A: S CHEMICAL POTENTIALS FOR STEP EDGE STRUCTURES

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).

### APPENDIX B: COMPLEX FORMATION ENERGIES

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_{m}S_{n}) > 0 of Cu and S atoms within the complex is given by^{20}

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_{m}S_{n} complex on terraces by extraction of Cu from the bulk satisfies E_{form}(Cu_{m}S_{n}) = δE(Cu_{m}S_{n}) − m *μ*_{Cu }− n δE(S), it follows that

Thus, E_{form}(Cu_{m}S_{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_{2} and Cu_{2}S_{3} via extraction from an extended Cu_{m}S_{m+1} chain for large m. From the definition of the S chemical potential, *μ*_{S}(Cu_{m}S_{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_{0} and select *μ*_{S} = *μ*_{S}(T) ≈ −1.91 eV and *μ*_{S}(chain) = −1.85 eV below. For extraction of CuS_{2} from the end of the chain, via the mechanism Cu_{m}S_{m+1} + S → Cu_{m−1}S_{m} + CuS_{2}, 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−1}S_{m} + CuS_{2} + Cu_{m′}S_{m′+1}, the formation satisfies

### APPENDIX C: CLUSTER EXPANSION FOR Cu_{m}S_{n} ENERGETICS

A perspective on the energetics of Cu_{m}S_{n} complexes follows from the application of a many-body cluster expansion commonly used in lattice-gas models for chemisorbed species.^{43} 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_{m}S_{n}) = −E_{form}(Cu_{m}S_{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_{2}) = 0.64 eV, and E_{bind}(Cu_{2}S_{3}) = 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_{3}S_{3} (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).