Graphene at Liquid Copper Catalysts: Atomic‐Scale Agreement of Experimental and First‐Principles Adsorption Height

Abstract Liquid metal catalysts have recently attracted attention for synthesizing high‐quality 2D materials facilitated via the catalysts’ perfectly smooth surface. However, the microscopic catalytic processes occurring at the surface are still largely unclear because liquid metals escape the accessibility of traditional experimental and computational surface science approaches. Hence, numerous controversies are found regarding different applications, with graphene (Gr) growth on liquid copper (Cu) as a prominent prototype. In this work, novel in situ and in silico techniques are employed to achieve an atomic‐level characterization of the graphene adsorption height above liquid Cu, reaching quantitative agreement within 0.1 Å between experiment and theory. The results are obtained via in situ synchrotron X‐ray reflectivity (XRR) measurements over wide‐range q‐vectors and large‐scale molecular dynamics simulations based on efficient machine‐learning (ML) potentials trained to first‐principles density functional theory (DFT) data. The computational insight is demonstrated to be robust against inherent DFT errors and reveals the nature of graphene binding to be highly comparable at liquid Cu and solid Cu(111). Transporting the predictive first‐principles quality via ML potentials to the scales required for liquid metal catalysis thus provides a powerful approach to reach microscopic understanding, analogous to the established computational approaches for catalysis at solid surfaces.

Experimental details. The graphene was grown in a dedicated CVD reactor designed for in situ characterization compatible with optical microscopy, Raman spectroscopy, and X-ray diffraction and reflectivity techniques [1]. A puddle of liquid Cu was used as a catalyst. For this, a few round pieces of 50-µm-thick Cu foils with a diameter of ~15 mm (99.9976% purity, purchased from Advent Research Materials, The United Kingdom) placed on a tungsten disk were melted and annealed above the melting point (1357 K) in Ar(91%)/H 2 (9%) atmosphere at a total pressure of 200 mbar (kept constant during the whole experiment) in the reactor prior to the Gr deposition. The Gr single layer was grown and characterized according to the protocol described in ref.
. [2]. The deposition was carried out in a 'pulse' mode. In this mode, the methane precursor(7 sccm of 2% CH 4 mixture with p CH4 =2.36 and p H2 = 7.45 mbar, respectively, in Ar for 20 s) was accumulated for some time before opening the valve and introducing the precursor into the reactor. After opening the CH 4 valve, the methane flow was stabilized at 7 sccm, equivalent to p CH4 /p H2 = 0.0127 (with p CH4 = 0.22 and p H2 = 17.19 mbar, respectively, at a total pressure of 200 mbar), until the layer was complete. The growth monitoring in real time was done through a quartz window on top of the reactor, where an optical microscope operated in radiation mode combined with a digital camera for image recording was installed.
The in situ XRR measurements were performed at the ID10 beamline of the ESRF synchrotron (Grenoble, France). To ensure the reactor transparency for the X-ray radiation with the beam energy of 22 keV, an aluminum side wall of the reactor wall was replaced with a copy made of thin beryllium. The vertical beam size was 15 µm, while horizontal size was 26 µm. A twodimensional MAXIPIX [3] detector was used. The XRR data from the curved liquid Cu surface were treated according to the method described recently in ref. [4] to eliminate the curvature's effect and extract the normalized reflectivity profiles. The resulting XRR curves were fitted using Refl1D software [5] to a slab model consisting of one (in case of bare Cu) or three (in case of Cu covered with Gr) layers: a Cu substrate, a void, and a C monolayer. The density of the liquid Cu, the C layer, and the thickness of the C layer were fixed to 7.99 g/cm 3 , [6] 5.36 g/cm 3 (derived from the graphite density [7], and 1.42 Å (diameter of the carbon atom [8]), respectively. The fitting parameters were the roughness of the Cu and Gr layers together with the width of the void layer. The roughness of the bare liquid Cu, of 1.1±0.1 Å, was obtained from fitting the reflectivity measured on it prior to Gr growth. This value was applied as a starting parameter for the fit of the Gr on the Cu curve. As a result, the roughness values of the Cu and Gr layer were 1.0±0.1 and 1.1±0.1 Å, respectively, and the void value was 1.5±0.1 Å. Note that both parameters (the roughness and the void width) are independent of each other.
Density Functional Theory Calculations. All the first-principles calculations were performed by the full-potential, all-electron DFT package FHI-aims [9] using light default settings for the basis set and integration grids, the PBE [10] or rSCAN [11] exchange-correlation functional, and a dipole correction. Dispersion corrections were included by either the Tkatchenko-Scheffler (TS) method [12] or via the MBD scheme [13]. K-point sampling was done with a density exceeding 60/Å k-points. The slab model for Cu(111) comprised 4 atomic layers with a (1x1) surface unitcell, which then hosted 2 C atoms of a slightly strained Gr sheet (Gr on Cu(111) in Table S1). The liquid Cu slab model has a comparable thickness. Electronic self-consistency converged energies to within 10 -6 eV, with the exception of large Gr-Cu(l) configurations computed with rSCAN, where a looser criterion of 10 -2 eV was used. Bulk fcc Cu and solid Cu(111) geometries were relaxed until residual forces fell below 1 meV/Å and 10 meV/Å, respectively.
Training of machine learning potentials The ML potential was trained based on the Moment Tensor Potential formalism [14] with hyperparameters of level=20 and a cut-off of 6 Å. The potential was trained against a dataset computed with different DFT approximations (see above). The training set includes multiple configurations corresponding to bulk and surface of the solid and liquid phases. Initial configurations were chosen from previous work [15] and then augmented by an iterative trial and error approach to improve the description of the interface in simulation cells tractable with DFT. The size of the final training and test set is 44 and 9, respectively. Multiple MTPs are trained by the same training set starting from different random seeds, and the final potential is selected based on both the training errors and the relaxed lattice parameters of fcc Cu and Gr. In addition, for every DFT approximation, the potential is retrained on a subset of the training set, which only contains bulk Cu and Cu surface (17 configurations). The effect of Cu-C training data on the Cu-Cu interactions in MTPs is eliminated in this way and the retrained potentials are only used for bulk Cu and Cu surface in Figure 7. The details of training and test configurations are summarized in Table S1 and all configurations are publicly accessible in Zenodo (https://doi.org/10.5281/zenodo.6993871).

Molecular Dynamics Simulations
The MD simulations are propagated for 1 ns with a time step of 1 fs. A Nose-Hoover thermostat [16] and a Parrinello-Rahman barostat [17] with damping parameters of 0.1 ps for temperature and 1 ps for pressure were employed to produce an NVT or NPT ensemble at various temperatures, respectively. To compute the adsorption height of Gr above liquid Cu, a cell composed by a liquid Cu slab (1589 atoms) and a Gr layer (364 atoms) is used. The simulation cell dimensions are adjusted to the 0K lattice constant for Gr. Since the liquid Cu slab can dynamically adjust its thickness and Gr has a small thermal expansion coefficient [18], this model system introduces minimal strain. For solid Cu(111), the simulation cell dimension matches the Cu lattice constant in order to avoid any strain in the metal slab. The Gr sheet is fitted into the simulation cell whereby an almost negligible strain of -0.1% is introduced in a cell that contains a solid Cu slab of 11664 atoms and a Gr layer of 3136 atoms. To evaluate the strain effect on the adsorption height, we build two further models for Gr-Cu(111), which contain 1680 atoms with +3.6% strain (1344 Cu atoms and 336 C atoms), and 11716 atoms with -0.5% strain (9216 Cu atoms and 2500 C atoms), respectively. All the solid Cu slabs have a thickness of 8 layers.
The structural data of these simulations were used to evaluate an averaged electron density profile along the surface normal (see Figure 3). To process the electron density profile, we smear a Gaussian function over the laterally-averaged atom density (see Figure S5). The area of the Gaussian corresponds to the number of electrons (29 for Cu and 6 for C), and the width σ is chosen similar to the experimentally determined roughness factors (σ =1.0 Å for Cu and σ =1.1 Å for C, see above). Specifically, we employed σ=0.9 Å for Cu and σ=1.1 Å for C, which led to the best match to the shape of the experimental profile. The small difference between experimental and theoretical width for Cu results from a fit via an error function in one case and a Gaussian in the other. We note that fitting the resulting electron density with an error function reproduces the choice of σ =1.0 Å for Cu in the experiment.
The inflection-point "gap" is determined as the distance between the average position of the C atoms and the inflection point of the Cu density profile (see Figure 2, 3 and S7). In order to determine the inflection point, we fit an error function to the electron density of Cu of which the inflection point and the half maximum coincide (see Figure S5b). This procedure is robust in comparison to numerically computing the second derivatives, which is prone to noise. We emphasize, that the σ for both the Gaussians and the error function only affects the electron density profile shape but not the determined Gr-Cu(l) gap if chosen reasonably large (see Figures S2 and S3). Furthermore, the averaged Gr-Cu "gap" is computed over MD trajectory segments of 0.1 ns over the last 0.9 ns and the uncertainties are derived from the standard deviation between these segments.
A 1 ns NPT MD trajectory of a liquid Cu bulk cell (864 atoms) at 1356 K (the experimental melting temperature of Cu) and one atmosphere of pressure provides the density of liquid Cu by every potential. For every potential, the NPT simulation repeats 3 times with different random seeds to get the mean value and the uncertainty of liquid density. The surface tension is evaluated from 1 ns MD simulations at 1370 K of a liquid Cu slab (1589 atoms) by this equation: where is the element of the pressure tensor and L is the length of the simulation cell along the z-axis. [19]. The way to compute the mean value and the uncertainty of the surface tension is just similar to the "gap". All the bond angles inside the Gr layer in the last 0.9 ns are collected to obtain the average bond angle to coordinating C atoms for every C atom and the signed delta angle is defined as = ( )(120 − ), where is the average bond angle of a C atom and is the z component of the displacement from the central atom to the centroid of its three neighbors. Here the sign of depends on whether the central atom is beneath or above its neighbors relative to the surface normal. The cell for the adsorption gap and energy of Gr on Cu(111) at 0K contains 4 Cu atoms and 2 C atoms and the C-C bond length in Gr is adjusted according to the relaxed fcc Cu (same as DFT calculations). The average adsorption energy at high temperature is obtained by the trajectory of the last 0.9 ns from the 1 ns simulation of the free-standing Gr, the liquid/solid slab and Gr on the slab. We used the atomic simulation environment (ASE) and LAMMPS for data handling and to perform simulations [20,21].

Comparison of the inflection-point and direct "gap"
To characterize the adsorption height of Gr on the solid Cu(111), we also compute the gap by a direct method. Here, the gap is determined as the distance between the (center of) the top surface atoms of the solid Cu slab and those of the Gr layer. Since the surface atoms are clearly defined for the solid slab, we can easily determine the average value and standard deviation of all the "direct" distances (see Figure S5) and sample them analog to the inflection-point "gap". This gap value corresponds to the experimentally determined distances for solid Cu [22,23] (and other metals). An inflection-point "gap" on the solid slab is in contrast not reported to our knowledge. Therefore, a straightforward comparison cannot be easily done based on literature data. On the other hand, such a direct "gap" for the liquid case is vague since a peak position is not sharp (see Figure S5d). Fortunately, through our simulations we can provide a comparison of the different "gaps" which is justified by the excellent agreement of both, the inflection-point gap for Gr on liquid Cu (see main text) and the direct gap for Gr on Cu(111). We determine this latter as 3.287(±0.0001) Å at 300K , which compares well to the experimental values of as 3.34 ± 0.06 Å [22] and 3.0 ± 0.2 Å [23] measured via ex situ total-reflection high-energy positron diffraction (TRHEPD) and atomic force microscopy (AFM) at room temperature. The direct value and the inflection-point determined gap in our simulations differ consistently by ≈ 1 Å at the employed values used for the electron density. This roughly stays constant with temperature and ML potential used. We therefore use this value (1 Å) to freely translate between the two "gap" definitions.