Finite-Temperature Thermodynamics of Cu(100) Oxidation: Missing-Row Reconstruction, Defect States, and Order-Disorder Transition from Nested Sampling

Abstract: Metal surfaces undergo structural, compositional, and morphological changes in response to their chemical environment. Tuning the surfaces' function and stability for a given application correspondingly necessitates an understanding of how this surface evolution couples to external conditions. Here, we demonstrate the feasibility of nested sampling simulations to obtain this coupling at first-principles predictive quality. By exploring the full configuration space, nested sampling estimates the partition function and gives direct access to desired thermodynamic ensemble averages at any temperature without prior knowledge. Computational feasibility is achieved through machine-learned interatomic potentials, an efficient GPU implementation of the sampling algorithm and bespoke sampling moves. Applied to the early oxidation of Cu(100), the approach successfully predicts the experimentally observed, complex $(2\sqrt{2}\times\sqrt{2})$R45$<sup>\circ$-O missing-row reconstruction. The full access to the partition function enables a detailed characterization of the temperature-dependent surface evolution, mapping the emergence of defect states and the order-disorder transition of the reconstructed surface.

• The paper extends nested sampling with a fine-tuned MACE potential, specialized Monte Carlo moves, and GPU acceleration to recover Cu(100) surface thermodynamics without assuming the reconstruction structure.
• Oxygen stabilizes the experimentally observed missing-row reconstruction, raises the simulated surface-melting temperature from 1060 K to 1400 K, and shifts Cu-only/Cu–O coexistence boundaries by up to 270 K through entropic effects.
• The study finds thermodynamically driven defect formation between 300 and 600 K and an order–disorder transition near 800 K, showing that the ordered reconstruction loses dominance well before surface melting.

Nested sampling has been established as a partition-function estimator for bulk materials, but its application to realistic surface systems has remained limited to model potentials. The work by Riccius, Reuter, Heenen, and Rogal extends nested sampling to a chemically complex surface problem—Cu(100) oxidation—and demonstrates that a single simulation can recover an experimentally observed reconstruction without prior structural assumptions while simultaneously resolving its full temperature-dependent evolution (2608.12787). The study combines nested sampling with a fine-tuned MACE machine-learned interatomic potential (MLIP), bespoke Monte Carlo (MC) moves for walker decorrelation, and a GPU-native implementation, and applies this machinery to the $(2\sqrt{2}\times\sqrt{2})\,R45^{\circ}$-O missing-row reconstruction (MRR) of Cu(100).

Methodological extensions

The central algorithmic challenge addressed is walker decorrelation under progressively tightening energy constraints $E < E_{\mathrm{limit}}$. Standard collective moves such as Galilean Monte Carlo (GMC) rely on small continuous displacements and fail once the energy ceiling falls below adsorption-site barriers: walkers become trapped in disconnected phase-space regions, biasing the ensemble. The authors augment GMC with lattice-based displacement moves—translations by integer or half-integer multiples of the surface lattice constant, the latter combined with a one-monolayer vertical shift to cross step edges—together with random displacements and identity swaps for the Cu–O system. A diagnostic on a single adatom over eight degenerate hollow sites shows that GMC alone produces strongly non-uniform site occupation, whereas the augmented move set restores uniformity.

Computational tractability is achieved by implementing the entire MC walk in PyTorch following the torch-sim framework, with batched neighbor lists and parallelized walks on a single GPU, yielding an approximately sixfold speedup relative to MPI-parallel ASE–MACE execution. Even so, cost remains substantial: roughly one week of wall time on an NVIDIA H100 for 24 mobile atoms and about one month for 40 mobile atoms.

The MACE potential was fine-tuned from a general Cu-surface-oxide model using active learning built around a genetic algorithm, adding 1200 Cu(100)-oxide configurations across six generations. It reaches test-set RMSEs of 1.7 meV/atom for energies and 45 meV/Å for forces. Simulations use a c$(4\times4)$ cell with three fixed substrate layers, a 15 Å mobile region bounded by a reflective wall, and 60 walkers per atom (1440 walkers for pure Cu, 2400 for Cu–O), with convergence enforced at 25 K via a negligible-probability criterion. Four independent runs per system show consistent heat capacities, supporting parameter adequacy.

Oxygen-induced reconstruction

Two compositions were studied at fixed $\theta_{\mathrm{Cu}} = 0.75$ ML: a Cu-only slab (24 mobile atoms) and a Cu–O slab ($\theta_{\mathrm{O}} = 0.5$ ML, 24 Cu + 16 O). All four independent runs converge to a single minimum-energy basin in each case. The pure-Cu global minimum is a heart-shaped anti-island with no missing-row character; the Cu–O global minimum is the ideal MRR observed experimentally. This contrast identifies the MRR as a purely oxygen-induced phenomenon, consistent with strain relief from oxygen intercalation: the Cu–Cu nearest-neighbor distance shifts from approximately 2.4 Å to 2.6 Å upon oxygen incorporation, and the second-neighbor RDF splits into three peaks reflecting distances parallel to, across, and within the Cu–O rows.

The heat capacity exhibits surface melting at 1060 K (pure Cu) and 1400 K (Cu–O), both below their respective bulk melting points as expected for surfaces, plus a pronounced low-temperature shoulder assigned to an order–disorder transition. The structure factor confirms long-range MRR periodicity only below roughly 800 K, well below melting, indicating that the reconstructed solid persists as an ordered phase over a limited window before disordering.

Defect states and the order–disorder transition

To classify the roughly 5.6 million sampled configurations, the authors employ canonical connectivity graphs with Cu–O edges below a 2.5 Å threshold—a descriptor robust against thermal vibrations—which resolves over 2.2 million topologically unique bonding patterns. Below 300 K the ideal MRR dominates exclusively. Between 300 and 600 K two defect states emerge: a bridging motif with two Cu atoms linking neighboring rows (partially formed perpendicular MR), reaching 22% population at 600 K, and a single-atom bridge intermediate reaching 8%. These motifs match the point-like and L-shaped patterns reported in STM studies after annealing at 500–600 K, where the predicted ideal-MRR population drops to 60%.

A key interpretive claim follows: these defects are thermodynamically competitive rather than kinetically trapped. This contradicts earlier attributions of such boundaries to step edges, domain mismatches, or domain-interconversion kinetics, and is consistent with computed diffusion and nucleation barriers below 1 eV that are surmountable above 400 K. Above typical annealing temperatures near 800 K—the heat-capacity shoulder—the entropy gain from defect proliferation outweighs the strain relief of the ordered MRR; the ideal MRR accounts for at most 20% of configurations, and beyond 1000 K all dominant structures fall below 5% probability, defining a disordered solid regime spanning roughly 600 K until surface melting at 1400 K. Some high-temperature configurations exhibit vibration-induced graph breakage, signaling incipient roughening.

Surface phase diagram

Absolute Helmholtz free energies from the partition function, combined with an O$_2$(g) reference chemical potential, yield a $(T, p_{\mathrm{O_2}})$ phase diagram analogous to ab initio thermodynamics but including full vibrational and configurational entropy. Relative to conventional zero-Kelvin energetics, the Cu-only/Cu–O coexistence line shifts upward by 30–270 K depending on pressure, quantifying entropic stabilization of the oxidized surface. The diagram further resolves intra-composition transitions: the ordered MRR disorders before oxygen desorption becomes favorable, and the molten Cu layer dominates at elevated temperature and low pressure. The authors note explicitly that the analysis is restricted to fixed compositions; intermediate phases reported previously (e.g., an ordered adlayer at $\theta_{\mathrm{Cu}} = 1.0$, $\theta_{\mathrm{O}} = 0.25$ ML) would modestly shrink the stability ranges found here, and grand-canonical extension remains open.

Limitations and open questions

Several limitations are conceded directly. The computational cost confines the study to two fixed compositions and a c$(4\times4)$ cell; larger cells and compositional variation require either faster potentials or grand-canonical nested sampling, which does not yet exist for surfaces. Walker counts (60 per atom) are slightly below those used for Lennard-Jones surfaces, and individual-run variability—including "extinction" events in which a secondary minimum basin loses all walkers—indicates finite resolution limits, partially mitigated by averaging over four runs. The origin of minor heat-capacity peaks near 1700 K (Cu) and 2350 K (Cu–O) remains unassigned, though structure-factor changes suggest subtle long-range-order modifications. Finally, the free-energy comparison treats $E < E_{\mathrm{limit}}$0 as a proxy for $E < E_{\mathrm{limit}}$1, valid only within the condensed-phase, low-pressure regime considered.

Conclusion

This work establishes nested sampling with MLIPs as a viable route to first-principles finite-temperature thermodynamics of non-trivial surface reconstructions. Applied to Cu(100) oxidation, it recovers the experimentally observed MRR as the global minimum without prior structural input, proves the reconstruction is oxygen-induced, and shows that its defect chemistry and order–disorder transition at approximately 800 K are thermodynamically driven rather than kinetic artifacts. The remaining obstacles are computational scale and ensemble extension, leaving open whether grand-canonical nested sampling can deliver composition-resolved surface phase diagrams at comparable accuracy.