Hydrogen Atom Scattering at the Al 2 O 3 (0001) Surface: A Combined Experimental and Theoretical Study

Investigating atom-surface interactions is the key to an in-depth understanding of chemical processes at interfaces, which are of central importance in many fields – from heterogeneous catalysis to corrosion. In this work, we present a joint experimental and theoretical effort to gain insights into the atomistic details of hydrogen atom scattering at the α -Al 2 O 3 (0001) surface. Surprisingly, this system has been hardly studied to date, although hydrogen atoms as well as α -Al 2 O 3 are omnipresent in catalysis as reactive species and support oxide, respectively. We address this system by performing hydrogen atom beam scattering experiments and molecular dynamics (MD) simulations based on a high-dimensional machine learning potential trained to density functional theory data. Using this combination of methods we are able to probe the properties of the multidimensional potential energy surface governing the scattering process. Specifically, we compare the angular distribution and the kinetic energy loss of the scattered atoms obtained in experiment with a large number of MD trajectories, which, moreover, allow to identify the underlying impact sites at the surface.


Introduction
Heterogeneously catalyzed reactions are of utmost economical importance, as a substantial part of the bulk products in chemical industry is made through reactions at solid catalyst surfaces 1,2 .Prominent examples are the Haber-Bosch process for ammonia synthesis 3 , the cracking of long-chain hydrocarbons in the petroleum industry 4 , and the oxidation of sulphur dioxide to sulphur trioxide finally yielding sulphuric acid 5,6 .The improvement of such catalytic systems crucially depends on a detailed understanding of the reaction mechanisms at the atomic scale that are governed by the potential energy surface (PES), which is a high-dimensional function providing the energy and forces for any given atomic configuration.Thus, the PES determines the dynamics of the reactants at the surface, the making and breaking of bonds and, finally, product formation.
Insights into the interactions between atoms or molecules and surfaces can be obtained from both, experimental as well as the-oretical studies.In particular atom and molecular beam experiments under well-defined conditions can provide invaluable details of processes at single-crystal surfaces [7][8][9] .For instance, information about the interaction potential and energy transfer between atoms or molecules and surfaces can be derived from inelastic scattering experiments.Moreover, using rare-gas atoms, surface phonons of a material can be studied 10,11 .Further, the role of rotational, vibrational and translational energy on surface reactivity is of central interest in molecule-surface scattering experiments 8,[12][13][14][15] .Finally, scattering of open-shell atoms like hydrogen or oxygen provides detailed information about chemical bond formation in surface reactions 9,16 .
Next to experiment, theoretical studies have become an indispensable tool in surface science and heterogeneous catalysis in the past two decades [17][18][19][20] .Still, in spite of the substantial increase in computational power, first-principles methods like density functional theory (DFT), which have been very successfully applied to the calculation of properties like adsorption energies and surface structures, remain too demanding for determining a large number of molecular dynamics (MD) trajectories when computing the forces on-the-fly for each new MD step.Thus, it is common practice to perform extended MD simulations for scattering at surfaces employing atomistic potentials, which -after fitting to accurate electronic structure data -can provide a direct relation between atomic positions and forces while they are much cheaper to evaluate.A wide range of methods has been proposed to represent atomistic potentials for surface scattering processes [21][22][23] , but constructing high-dimensional PESs of first-principles quality taking all the surface degrees of freedom explicitly into account has remained a substantial challenge.
In recent years, the introduction of machine learning potentials (MLP) [24][25][26][27][28][29] has led to a paradigm change in the development of PESs, since flexible machine learning algorithms allow to combine the accuracy of electronic structure methods with the efficiency of simple empirical potentials.In fact, molecule-surface interactions have been among the early applications of MLPs [30][31][32][33][34] , which in the first years still employed a frozen-surface approximation to restrict the complexity of the PES.Starting with the introduction of high-dimensional neural network potentials (HDNNP) in 2007 35,36 , the construction of MLPs for high-dimensional condensed systems has become possible.Nowadays, many different types of MLPs are available [37][38][39][40][41][42] , which allow to explicitly take all degrees of freedom into account even for systems containing thousands of atoms.This methodical progress now enables to study atom-and molecule-surface interactions with full dimensionality including a mobile surface, and numerous applications have been reported to date [43][44][45][46][47] .
In spite of these advances, so far most studies have addressed the surfaces of rather simple materials like metals 44,[48][49][50] , solid xenon 51 or graphene 43,52 , i.e., pure elemental systems, while scattering at binary compounds like oxides or thin surface oxides 53 has been rarely studied to date.Only recently, a first Hatom beam experiment has been reported for the interaction with α-Al 2 O 3 54 .Its surface termination is well-studied under high vacuum conditions and does not show any surface reconstruction.The material is also often applied, either as catalyst 55 or as supporting material 56 , increasing the need for a deeper understanding of processes at its interfaces.
Here, we report a combined theoretical and experimental study addressing the interaction of H-atoms with the α-Al 2 O 3 (0001) surface with the aim to unravel the atomistic details of the scattering process.For this purpose we perform H-atom beam experiments employing four different incident conditions to probe the scattering angle distribution as well as the kinetic energy loss of the atoms leaving the surface.Unlike other systems, such as those involving vibrations or steric effects of scattering molecules, in the present work on an atomic beam only the translational energy is involved.We study the effect of this energy for two different incidence angles.For a deeper analysis, in parallel, we carry out large-scale molecular dynamics simulations based on energies and forces obtained from a HDNNP, which after training provides a first-principles quality description of the scattering process.This combination of methods is used for a detailed analysis of the quality of the PES for different initial beam kinetic energies and incidence angles, which allows us to assess the reliability of the employed exchange correlation functional for scattering at different α-Al 2 O 3 surface sites.

Experimental methods
The hydrogen atom scattering instrument employed in this work is described in detail in Refs. 9,57.It utilizes a mono-energetic H-atom beam formed by photolysis of a supersonic jet of hydrogen iodide molecules using an excimer laser operating with KrF (248 nm).The generated H-atoms pass through a skimmer and two apertures separating two differential pumping regions and hit the sample in the ultrahigh vacuum (UHV) scattering chamber.The sample is mounted on a 6-axis manipulator that allows the variation of the polar and azimuthal incidence angles.Scattered H-atoms are detected by Rydberg atom tagging timeof-flight (TOF) 58 , where two laser pulses excite the H-atoms to a long-lived (n = 34) Rydberg state.The neutral Rydberg atoms travel 250 mm, pass a grounded mesh and are field-ionized and detected with a multichannel plate (MCP) detector.A multichannel scaler records the TOF distributions that are converted to kinetic energy distributions.The detector is rotatable in the plane containing the H-atom beam and the surface normal, providing TOF spectra at various scattering angles.The alumina sample was cleaned by annealing for several hours in an oxygen atmosphere (10 × 10 −6 mbar) at 600 • C. The cleanliness and structure of the surfaces are monitored by Auger electron spectroscopy (AES) and low-energy electron diffraction (LEED).Based on AES and LEED we conclude that the surface has a (1×1) structure with Al termination, Al-O 3 -Al trilayers [59][60][61][62] and a step density of about 2-4%.Moreover, it is important to note that depending on the stacking of the trilayers there are two possible terminations related by a mirror operation resulting in two possible incident azimuthal angles φ i with respect to the [10 10] surface direction of either 0 • or 180 • 63 that cannot be distinguished with the methods available in experiment, which needs to be considered in the simulations.Further details about the sample characteristics are given in the SI.The inelastic scattering of the H-atoms was performed for initial kinetic energies of 0.99 eV and 1.92 eV and incidence polar angles of θ i = 40 • and θ i = 55 • resulting in four different sets of scattering conditions.

High-dimensional neural network potentials
The full-dimensional potential energy surface of the hydrogen atom scattering process at α-Al 2 O 3 is represented by a highdimensional neural network potential.HDNNPs are a frequently used type of machine learning potential introduced by Behler and Parrinello in 2007 35 .In second-generation HDNNPs as employed here 64 , the total energy of the system is constructed as a sum of atomic energies, that depend on the local atomic environments up to a cutoff radius R c , which has to be chosen large enough to include all relevant atomic interactions.The positions of all neighboring atoms inside the resulting cutoff spheres are described by vectors of atom-centered symmetry functions (ACSF) 65 , which are invariant with respect to rotation, translation and permutation of the system and in combination with Eq. 1 allow to explicitly include all atomic degrees of freedom.The atomic energies are then computed as outputs of element-specific atomic neural networks as a function of the respective input ACSF vectors.Since only the Cartesian coordinates of all atoms in the system are required for the calculation of the ACSFs, HDNNPs are able to describe the making and breaking of bonds, which is an essential requirement for studying scattering at surfaces.
The weight parameters of the atomic neural networks are obtained in an iterative training process making use of energies and forces obtained from electronic structure calculations 66,67 .Since atomic energies are not quantum mechanical observables, total energies are used as target properties, and the partitioning into atomic energies is done automatically during the training process.The forces, which are needed for training as well as for running molecular dynamics simulations, can be computed as analytic derivatives of the HDNNP energy.More details about HDNNPs, the underlying methodology and typical applications can be found in several reviews 36,64,66,68 .
3 Computational Details

Density functional theory calculations
The DFT reference energies and forces for training the HDNNP have been calculated using the Fritz-Haber-Institute ab initio molecular simulations (FHI-aims) code (version 171221_1) 69 applying the "light" settings for the integration grid and for the basis sets, which consist of numerical atomic orbitals.Spin-polarized calculations were performed for structures containing hydrogen atoms, while α-Al 2 O 3 bulk and slab structures without hydrogen atoms have been treated as closed-shell systems in a numerically consistent way.Γ-centred k-point grids have been used for all periodic systems.For the bulk α-Al 2 O 3 structures containing 40 atoms a 6×6×2 k-point grid was used, while for the (2×2) slab supercells a 2×2×1 k-point grid has been employed.Using these k-point grids, the formation energy of bulk α-Al 2 O 3 and the adsorption energies of hydrogen atoms at different surface sites are converged to about 0.5 meV.The convergence criteria for the electronic self-consistency of the single point calculations have been set to 10 −6 eV for total potential energies and 10 −4 eV Å −1 for the forces.
Several exchange correlation functionals have been tested in combination with Tkatchenko-Scheffler dispersion corrections 70 , including the GGA functionals PBE 71 and RPBE 72 as well as the PBE0 73 hybrid functional.Figure S3 in the SI shows two onedimensional cuts through the potential energy surface for a hydrogen atom approaching the surface on top of an aluminium and an oxygen surface site, respectively.For both sites the the observed interaction energies at the potential minimum are within a range of 0.1 eV for all functionals, which is the expected uncertainty with respect to the description of exchange and correlation.This is about one order of magnitude smaller than the lower initial kinetic energy of 0.99 eV used in experiment, even if a potential energy amount of 0.1 eV would be converted entirely to H-atom kinetic energy.At much lower experimental kinetic energies, however, the accuracy of the employed functional would become highly relevant.At both surface sites the PBE functional is more attractive than RPBE with the minimum region showing no clear trend for a better agreement with the hybrid functional.Hence, in the present work, which strongly depends on the de-scription of the repulsive walls due to the high kinetic energies of the H-atoms, we will employ the RPBE functional that is in reasonable agreement with PBE0 for smaller, i.e., repulsive, atomsurface distances.Performing a very large number of PBE0 calculations would be very demanding for the construction of the training set while on the other hand using this hybrid functional does not yield a qualitatively different description of the PES.

Construction of the reference set
The reference data set contains bulk geometries and slab structures of the α-Al 2 O 3 (0001) surface, both with and without an H-atom at various positions above the slab.The bulk structures contain one unit cell of α-Al 2 O 3 , while for the slabs (2×2) supercells of the α-Al 2 O 3 (0001) surface unit cell have been constructed with the vacuum aligned in z direction to avoid lateral interactions between the periodic images of the H-atoms.The slabs consist of 8 Al-O 3 -Al trilayers as shown in Figure 1.The surfaces are terminated by one Al-atom per unit cell on top of an oxygen layer, which is called the "half metal" termination and corresponds to the surface structure present in experiment [60][61][62] .The atoms of the bottom four trilayers have been frozen at their bulk equilibrium positions, while the atoms in the top four layers have been relaxed and are also mobile in the MD simulations, which ensures that the frozen surface atoms in the deep layers are outside the cutoff spheres of the ACSFs describing the hydrogen atom environment.Utilizing an analysis of the atomic interactions based on the Hessian matrix of the slab 75 , we found that the remaining interactions between the deep frozen layers and the H-atom are negligibly small even for H-atom positions close to the surface (see supporting information).The vacuum between the slabs was set to 21 Å.At a distance from the surface of about 6.3 Å corresponding to the ACSF cutoff only very small forces of roughly 0.003 eV Å −1 act on the hydrogen atoms such that atomic interactions at larger separations can be neglected to a good approximation.The functional form of the cutoff function 65 ensures a smooth decay of the forces to zero at the cutoff radius in the HDNNP.The size of the vacuum prevents any notable interactions between the H-atom and the opposite side of the slab in the periodic DFT calculations and the HDNNP even for the largest H-atom-surface distance of 6.8 Å we studied.
The reference data points were generated combining various approaches.Initially, 133 bulk and 200 clean slab structures were generated with random atomic displacements of up to 0.2 Å.Initial HDNNPs based on these structures were used in MD simulations in the canonical ensemble at increasing temperatures up to 400 K to identify important missing surface structures through active learning [76][77][78] .Structures showing an energy or force variance larger than 2.8 meV atom −1 or 0.5 eV a 0 −1 (with a 0 being the Bohr radius) for different preliminary HDNNPs, respectively, were added to the data set.
Hydrogen atom positions above the surface were initially included using an equidistant grid with spacing of 0.5 Å in the x and y Cartesian directions, starting at the vertical position of the topmost oxygen atom layer and up to 6.75 Å above this layer within the symmetry-unique part of the unit cell (s.Fig. 1) employing the fixed structure of the clean relaxed surface.In the z direction the distances between the H-atom positions are increasing with distance from the surfaces due to the reduced variance in the atomic interactions at larger distances.Up to a distance of 1.0 Å from the topmost oxygen layer a vertical grid spacing of 0.125 Å has been employed, followed by a vertical spacing of 0.25 Å for distances up to the limit of 6.75 Å.From this grid, all strongly repulsive structures with a distance between the H-atom and any surface atom below 0.5 Å have been discarded.In the next step, active learning has been performed employing molecular dynamics simulations of the scattering process using the experimental incidence conditions including a mobile surface to sample the surface degrees of freedom.Additional structures obtained in these trajectories were additionally selected by farthest point sampling 79 as discussed in the SI to improve the representation of rarely visited geometries.

Construction of the high-dimensional neural network potential
For the construction of the HDNNP the RuNNer code (version from August 22, 2019) was used 36,66 .The atomic neural networks consist of two hidden layers containing 19 neurons each with n α G input neurons for atoms of element α corresponding to the respective number of ACSFs, and one output node providing the atomic energy.The ACSF cutoff radius is R c = 12 a 0 (6.35 Å), which is sufficient to include a major part of the atomic interactions as discussed in Section 3.2.A list of the employed atomcentered symmetry functions can be found in the supplementary information.Since there is only one H-atom per structure, ACSFs for the description of H-H interactions have not been included.For the construction of the HDNNP the atomic energies were removed from the total energies before training such that numerically more favorable binding energies are learned.The HDNNP was then trained to reproduce these binding energies and the DFT atomic force components of the reference structures.The reference data set was randomly split into a training set consisting of about 90 % of the structures used for adjusting the weight parameters and a testing set of the remaining 10 % of the structures for assessing the quality of the HDNNP for unknown structures.The parameters of the Kalman filter 80 were set to λ = 0.98000 and ν = 0.99870.Overall, the structures cover an energy range of 0.446 eV atom −1 .Atoms with forces up to 13.6 eVÅ

Molecular dynamics simulations
To run MD simulations the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (version from 16 th March 2018) 81 including the n2p2 extension for HDNNPs 82 was used.MD simulations for active learning were run in the canonical NV T (bulk and clean surfaces) and in the microcanonical NV E (trajectories for H-atom scattering) ensemble with a time step of δt = 0.5 fs for simulations without a hydrogen atom and δt = 0.25 fs for simulations including a hydrogen atom.The velocity Verlet algorithm was used as integrator 83 .To simulate the NV T ensemble a Nosé-Hoover thermostat 84 with a temperature damping parameter of 0.01 ps was used.
For the MD trajectories of H-atom scattering, i.e., slab structures at room temperature, atomic velocities of the surface atoms were taken from NV T trajectories at 300 K in intervals of 0.1 ps.Then, the H-atom was placed at a random lateral position 6.3 Å above the surface with velocity vectors corresponding to the experimental kinetic energy, polar and azimuthal angle.The trajectories were run in the NV E ensemble and terminated after 0.5 ps or before, when the H-atom surface distance exceeded a value of 7.8 Å. Trajectories were considered in the analysis of the kinetic energy loss and the angular distributions if the H-atom was at least at a distance of 7.3 Å from the slab at the end of the trajectory.The resulting kinetic energies and scattering angles of the H-atoms were calculated using the velocity vector of the last time step of the trajectory.The detector counting the H-atom flux is defined as a vector with the specified polar angle θ s towards the surface normal and in plane with the incident azimuthal angle φ i to mimic the position of the detector in the experiment.The experimental detection angle γ, which is the angle between the velocity vector of the H-atom and the detector, is about 1.5 • , while in the MD simulations larger angles up to 5 • have been included to improve statistics, which has very little effect on the obtained results (for a discussion see SI).

Phonon Calculations
Phonon band structures were calculated with the phonopy code (version 2.17.1) 85,86 using both, force constants from FHI-aims RPBE DFT calculations as well as force constants calculated using the HDNNP.The force constants were determined using a geometry optimized primitive α-Al 2 O 3 cell containing 10 atoms.The geometry optimization using the HDNNP was carried out using LAMMPS 81 with the n2p2 82 extension.To determine the correct symmetry, the symmetry precision threshold was set to 5 • 10 −4 .The interface between phonopy and RuNNer was created using a Python script.respect to the reference method are significantly smaller than the uncertainty with respect to the choice of the exchange correlation functional.More detailed information about the energy and force errors can be found in the SI.

Results and Discussion
As monitoring the RMSEs of the energies and forces alone is insufficient to judge on the quality of a potential, in a first validation step we calculated the lattice constants and bulk modulus B 0 for bulk α-Al 2 O 3 .As can be seen in Table 1, the HDNNP predictions are in excellent agreement with the underlying RPBE-DFT data, while compared to experiment the RPBE functional slightly overestimates the lattice constants and underestimates the bulk modulus.In a next step, we have computed the phonon band structure of bulk α-Al 2 O 3 using the phonopy program 85,86 .The phonon spectrum obtained from the HDNNP in Figure 2 shows excellent agreement with DFT in particular at low frequencies.Overall, the phonon band structure, which is important for energy dissipation after the scattering process, is well represented, and we conclude that the HDNNP provides a reliable description of bulk α-Al In order to ascertain the quality of the structural description of the clean surface, we compared the relaxed interlayer distances of the clean surface optimized with RPBE DFT and with the HDNNP.When optimizing the surface with DFT, the topmost three interlayer distances relax from the initial bulk separations of (Al-O 3 ), 1.69 Å (O 3 -Al) and 0.52 Å (Al-Al).Consequently, there is a substantial inwards relaxation in particular of the top Al layer, and in the HDNNP relaxation this Al-O 3 layer distance is about 0.15 Å with an error of only about 0.016 Å.For all other interlayer distances the differences between the HDNNP and RPBE results are below 0.01 Å.
As an initial assessment of the quality of the hydrogen atomsurface interactions, we have analyzed the 3D-PES of a H-atom at different positions above the surface, which, as described in Section 3.2, also formed the starting point of the reference data generation beyond this grid.Figure 3 shows several 2D cuts through this PES at different H distances from the surface, which is frozen in its relaxed geometry for this purpose.Overall, the RPBE and the HDNNP PESs displayed in Figures 3a and 3b are very similar, as can also be seen in the energy difference plot shown in panel (c).The total energy differences are found to be typically between −0.01 eV and 0.01 eV, which corresponds to an error of 0.06 meV atom −1 .The largest difference is about −0.04 eV in the cuts at 2.0 Å and 2.25 Å, i.e., 0.25 meV atom −1 , which is still well below the RMSE of the HDNNP.We conclude that the differences are very small.However, the 3D-PES only probes the interaction between the H-atom and the frozen surface, which is insufficient to check the reliability in MD simulations of the scattering process with fully mobile surface atoms.
To validate the performance of the HDNNP for full-dimensional trajectories, ab initio MD simulations have been carried out and compared to the HDNNP.The starting structures for the ab inito MD simulations and the atomic velocities for the slab atoms were taken from NV T simulations at 300 K of the clean surface as determined using the HDNNP.The magnitude and direction of the velocities of the H-atoms correspond to the experimental conditions, and we have chosen such conditions that are estimated to result in experimentally detectable scattering events based on preliminary HDNNP trajectories.In total, ten different ab initio MD trajectories have been computed, and the results are shown in Table 2. Since even smallest numerical differences result in diverging MD trajectories, a direct comparison of independently computed ab initio and HDNNP trajectories using the same initial conditions is not possible.Thus, for assessing the reliability of the HDNNP, we have recomputed the structures visited in the ab initio trajectories by single point HDNNP calculations and computed the RMSEs of the energies and the z-components of the forces acting on the H-atom (s.Table 2).Overall, we find a very good agreement between DFT and the HDNNP for all trajectories with deviations similar to or below the RMSEs of the training and test data sets.
Figure 4 shows a detailed analysis of the scattering process for the example of trajectory MD4 with an incident kinetic energy of 1.92 eV and an incident polar angle of 40 • .The trajectory was calculated using RPBE DFT and afterwards the structures were recalculated using the HDNNP.A comparison of the DFT and HDNNP potential energies of the system along the trajectory is shown in Figure 4a, while Figure 4b displays a comparison of the DFT and HDNNP z-components of the forces acting on the H-atom.Both, energies and forces, are very well represented by the HDNNP along the entire trajectory, which is shown in a top and a side Table 2 Initial conditions of the ab initio MD simulations used to validate the HDNNP, including the incident polar angle θ i , the incident azimuth angle φ i and incident kinetic energy E kin,i , as well as the polar angle θ s and the kinetic energy of the scattered atoms E kin,s .Moreover, the RMSE values of the total energy and of the z-component of the H-atom force vector have been computed by comparing the ab initio MD trajectory with subsequent HDNNP single point calculations along the trajectories to quantify the deviations between DFT and the HDNNP.The surface has been equilibrated at 300 K.  and d).We note that the energy deviations are small in comparison to the initial kinetic energy of the collision of 1.92 eV and thus do not significantly affect the scattering process.
A similar agreement has been found also for the other trajectories in Table 2.

Comparison of experimental scattering and simulations
Having validated the accurate description of the DFT-PES by the HDNNP, we now perform large-scale MD simulations and compare the scattering properties of the α-Al 2 O 3 (0001) surface with experiment using all four scattering conditions as described in Section 2.1.Due to the close agreement between DFT and the HDNNP, a comparison between experiment and simulation will allow to assess the performance of the employed DFT description of the scattering process.Concerning the experimental side, in view of the rather high complexity of the α-Al 2 O 3 (0001) surface in comparison to previous atomic beam studies at surfaces of pure elements, we took great care when characterizing the surface prior to the scattering experiments in detail (see SI).Still, we cannot exclude that differences between experiment and simulation may be related to experimental shortcomings, and a comparison of experimental and theoretical results needs to be done with care.
In the experiment, incidence polar angle, incidence kinetic energy and surface temperature are well-controlled.Also the incident azimuthal angle is experimentally defined, but here we have to consider that the surface can be terminated in two ways as discussed in Section 2.1.However, both represent very different surface morphologies for a tilted incidence angle of the beam and thus produce significantly different scattering spectra.Due to the lack of more detailed information, for the comparison between experiment and MD we have computed trajectories of both terminations assuming a 1:1 ratio in experiment, which in fact also provides the best overall agreement.A further discussion of this aspect can be found in the SI.
Figure 5 shows the experimentally obtained kinetic energy loss spectra for the incidence polar angle of 40 • and both initial kinetic energies E kin =0.99 eV and 1.92 eV along with the respective data obtained in the MD simulations.The corresponding figure for the incidence polar angle of 55 • is given in the SI.The distributions have been normalized to a maximum peak height of one.For all incidence conditions a sticking probability of at most 5 % has been found (s.Table 3), which therefore does not significantly affect our results.Overall, it can be seen that the simulations have a slight tendency to overestimate the energy loss of the scattered H-atoms with respect to experiment, achieving better agreement for the larger incidence energy.
Figure 6 compares energy integrated angular distributions of the scattered H-atoms obtained from experiment and simulation.The distributions have been normalized to a maximum of one and the theoretical distributions are again combinations of the two possible surface terminations φ i =0 • and 180 • .For all incident conditions the position of the maximum compares well to the experiment and there is very good agreement between experiment and theory for large scattering angles.However, at smaller scattering angles experiment and theory exhibit some deviations, as in experiment scattering angles close to the surface normal are more frequently found compared to the MD simulations.

Discussion
Although overall a reasonable agreement has been obtained, there are some deviations between the experimental data and their counterparts determined by molecular dynamics as in the  simulations there is a larger kinetic energy loss (Fig. 5) in particular for the lower incidence kinetic energy as well as a narrower angular distribution of the scattered atoms (Fig. 6) compared to the atomic beam results.However, in previous work, good agreement between experiment and simulation was achieved for welldefined, single-element surfaces, i.e., solid Xe 51 , graphene 43 and several fcc metal surfaces 49 .The more pronounced discrepancies in the present study may have several reasons.First, it has to be ensured that the surface structure used in the simulations is a good representation of the surface that is present in experiment.In case of the α-Al 2 O 3 (0001) surface, controlling the surface structure is more challenging, as, e.g., there are two possible terminations, which are geometrically inequivalent for the scattering atoms and different ratios of both terminations is essentially unknown.The effect of this ratio on the simulation outcome is discussed in the SI.Moreover, the effect of steps and other surface imperfections is not included in the simulations.
Concerning the simulations, the accuracy of the reference electronic structure method poses a limit for the quality of the HDNNP that can be achieved.In the present work a GGA functional has been employed as our initial investigations have shown that for the investigated geometries there are no fundamental differences between the quality of a GGA description and the PBE0 functional in view of the high atomic kinetic energies.However, due to the substantial computational costs of the hybrid functional it has not been possible to construct a full-dimensional HDNNP based on PBE0, and consequently parts of the PES exhibiting larger deviations might have been missed.
We will start discussing the results of the simulations before we address the deviations between experiment and theory in more detail.A first step to gain some basic information about the scattering mechanism is to compare the average energy loss with the predictions of the simple Baule model.The Baule limit describes the maximum energy loss that can be expected in a single-atom surface collision.It is based on a zero impact parameter collision of the impinging H-atom with a resting surface atom.An energy loss significantly larger than the Baule limit is a clear indication for a complex scattering mechanism, for example multiple collisions or energy transfer to several atoms in a single collision.Table 4 compares the Baule limit for the oxygen and the aluminum site to the average energy losses observed in experiment and theory.It can clearly be seen that the simulation, especially for an incidence kinetic energy of 0.99 eV, yields a larger relative energy loss than the simple Baule limit, while in experiment a loss below the Baule limit is observed.This observation indicates that there might be a complex scattering mechanism in the simulations that is not necessarily expected based on the experimental data.
To understand the reason for the large energy loss observed in simulation we took a closer look at the scattering trajectories leading to high energy losses and analysed them in detail.First, we extracted the geometry of closest approach of the H-atom to the surface for each trajectory to allocate the surface impact sites and to correlated them with the energy loss and scattering angle.In Figure 7 the closest point of approach for one incident condition, θ i = 40 • and E kin = 1.92 eV, is plotted.The black lines highlight a surface unit cell of α-Al 2 O 3 .Although the surface is mobile at the surface temperature of 300 K, the estimated location is accurate enough to analyse the origin of the bounce and the trajectories show a clear trend.In trajectories leading to small scattering angles the H-atom is mainly scattering from the topmost Al atoms, at medium angles the scattering predominately takes place at the O atoms and hollow sites between them, while at the largest scattering angles the scattering is again mainly occurring at the Al atoms.This scattering behavior can be explained by rather simple geometric considerations.In Figure 8 a schematic representation of elastic scattering processes at different surface sites is shown for an incidence angle of θ i = 40 • .The top-layer Al atoms effectively shield the O layer preventing large and small scattering angles for the case of scattering at oxygen atoms, while collisions with the Al atoms can result in very low and large scattering angles.The rather flat oxygen layer thus only allows for scattering at angles close to the medium-sized incidence angle.These geometric conditions also apply to the other investigated incident conditions.Best agreement between simulation and experiment for E kin = 1.92 eV is achieved for small and large scattering angels for which scattering from Al dominates.For intermediate scattering angles, where scattering from O dominates, the simulation results deviate more strongly from experiment.This might be a hint for a less accurate description of the scattering process at the oxygen sites.Since the comparison of the HDNNP trajectories with ab initio MD does show a similar good agreement for oxygen and aluminium sites, the most likely explanation thus is a less accurate description of the hydrogen-oxygen interaction by the RPBE functional.
Having identified that trajectories initially scattered at oxygen atoms tend to show a too large energy loss, we took a closer look at representative trajectories to get a better understanding of its origin.Specifically, we analyzed the atomic kinetic energies and the changes of the positions of the surface atoms in selected trajectories by comparing them to trajectories with the same initial surface atom velocities but without the presence of a H-atom.In this way, we are able to see how all surface atoms are influenced by the collision with the H-atom.In this analysis we observed that if the H-atom hits an oxygen atom on the surface, the kinetic energy is transferred not only to one oxygen atom but is rapidly distributed over several neighboring surface atoms.Due to the involvement of multiple indirect collision partners during the collision with the O-atoms, the energy loss is increased compared to an idealized single-atom collision.Figure S15 in the SI illustrates how the kinetic energy of the H-atom is efficiently transferred to multiple surface atoms.A similar mechanism has been observed for H-atom scattering from graphene 52 .If, however, the H-atom is initially scattered at an Al-atom, the energy is mostly trans- ferred to this Al with an overall much smaller energy loss in a more elastic process.This efficient energy transfer is also seen in ab-initio MD simulations supporting that this effect is not a consequence of inaccuracies in the HDNNP.
One possible source of error contributing to deviations between experiment and simulation are surface imperfections.Although we took great care in experiment to prepare a well-defined surface, some defects are unavoidably present.From AFM measurements we know that the step density of the surface is about 2-4%.Furthermore, it is known that the α-Al 2 O 3 (0001) surface loses oxygen in an UHV environment.A significant O loss would lead to surface reconstruction and can be excluded, since we observed the (1×1) LEED pattern at all times.Still, oxygen vacancies will be present at the surface.Furthermore, the presence of hydrogen on the surface can not be excluded 60 and may affect the energy loss of the scattered hydrogen atoms.The deviation in the angular distributions is a clear indication that the surface is not perfectly flat.Steps and defects increase scattering normal to the surface, exactly what we observe (s. Figure 6).However, our finding that the energy loss in experiment is smaller than in simulations can not easily be assigned to surface defects, as surface defects commonly lead to an increased energy loss, for example due to a higher probability of multi-bounce scattering.However, since a complex scattering mechanism involving many atoms is causing the larger energy loss at the oxygen sites in the simulations, it cannot be excluded that surface defects could change or even prohibit such a process.
The second point to be addressed is the general accuracy of the RPBE exchange correlation functional, which determines the interaction strength in between the surface atoms and between the surface atoms and the H atom.The best agreement between the experimental and theoretical kinetic energy distributions is obtained at low scattering angles with high kinetic energies.The largest difference in mean kinetic energy loss is present at large scattering angles and low kinetic incident energies, especially the scattering at E kin,i = 0.99 eV and θ s = 65 • .Both cases represent extreme cases of H-atom-surface interaction.When the H-atom scatters at high kinetic energies and low scattering angles, the effective interaction time is minimal.However, when the scattering angle is large and the kinetic energy is low, the interaction time with the surface is significantly longer.As the kinetic energy distributions for high kinetic energies are in better agreement with the experiment compared to the distributions for low kinetic energies, apparently inaccuracies of the reference electronic structure method tend to accumulate with interaction time.Moreover, the bulk modulus of α-Al 2 O 3 obtained in the DFT calculations is about 10% smaller than the experimental value (Table 1).This softer description of α-Al 2 O 3 by the RBPE functional may lead to an overestimation of the kinetic energy loss.

Summary
We have studied the interaction of hydrogen atoms with the α-Al 2 O 3 (0001) surface using a combination of large-scale molecular dynamics simulations based on a high-dimensional neural network potential trained to DFT data employing the RPBE functional and highly accurate H-atom beam scattering experiments under UHV conditions.Two different initial kinetic energies and two different incidence polar angles have been investigated.Best agreement between experiment and theory is found for large initial kinetic energies and very low and high scattering angles, which we find to arise from scattering at top-layer aluminium atoms.Scattering at lower initial kinetic energies results in a larger loss of kinetic energy in the MD trajectories compared to experiment, and in general also scattering at oxygen sites seems to result in larger discrepancies between experiment and theory.As a careful validation of the multidimensional PES representation by the HDNNP showed a very close agreement with the underlying DFT reference calculations, we consider limitations in the accuracy of the employed RPBE functional as the most likely explanation of the observed deviations.Although, in view of the high kinetic energies, the differences between different exchange correlation functionals are small, more complex scattering mechanisms at the oxygen sites might increase the sensitivity Notes and references

Fig. 1
Fig. 1 Side (a) and top view (b) of the relaxed (2×2) α-Al 2 O 3 (0001) supercell slab.Al atoms are shown in grey, O atoms are displayed in red.The white atoms represent the regular H-atom grid employed for the initial screening of the potential-energy surface using a frozen surface.The orange lines in (b) represent the boundaries of the symmetry-unique wedges, and two examples for the three-fold rotational symmetry axes aligned perpendicular to the surface are highlighted.The figures have been created using Ovito version 3.8.4 74.

− 1
have been included in the data set for monitoring the training errors, but only forces up to 3 eVÅ −1 have been used to update the weights in the training process.

4. 1
Validation of the HDNNP The final data set obtained after several active learning cycles consists of 15,812 structures, which include 808 bulk configurations and 15,004 slab geometries (12,704 with H-atom and 2,300 clean surface structures).The final HDNNP has a root mean square error (RMSE) of 0.746 meV atom −1 for the energy and 0.103 eV Å −1 for the atomic force components for the test data set, which is in the typical order of magnitude of state-of-the-art MLPs.For the training set the RMSE values are 0.257 meV atom −1 and 0.111 eV Å −1 , respectively.Therefore, the HDNNP errors with

Fig. 2
Fig. 2 Phonon band structure of bulk α-Al 2 O 3 computed for the HDNNP and DFT employing the RPBE functional using the phonopy program 85,86 .

Fig. 3
Fig. 3 2D cuts at different distances through the PES of the H-atom above a frozen α-Al 2 O 3 (0001) surface.Panel (a) depicts the (2×2) supercell of the RPBE PES while (b) shows the HDNNP PES.Given are relative potential energies ∆E of the entire system with respect to the energy of a H-atom at infinite separation from the surface.The underlying data points correspond to the regular grid of H-atom positions shown in Fig. 1.Panel (c) shows the difference ∆∆E between RPBE and the HDNNP.Note that in contrast to the energy RMSEs these energies are not normalized per atom in the system.

Fig. 4
Fig. 4 Analysis of ab initio MD trajectory MD4 (see Table2).Panel (a) shows the DFT and HDNNP potential energies of the structures ∆E slab+H with respect to the relaxed slab and the H-atom at infinite distance from the surface vs. the time t of the trajectory.The trajectory was calculated using RPBE DFT and afterwards the structures were recalculated using the HDNNP showing very close agreement.Panel (b) shows the z-component of the force vector acting on the H-atom.Panels (c) and (d) show the path (blue) of the H-atom (white) along the surface (top and side views).The surface structure corresponds to the first frame of the simulation, the highlighted H-atom denotes the starting position above the surface.

Fig. 5
Fig. 5 Kinetic energy loss distributions for H-atom scattering at a) E kin,i =0.99 eV and θ i = 40 • , b) E kin,i =1.92 eV and θ i = 40 • .The black line shows the experimentally measured distributions, while the red line shows the theoretical distributions.The detected scattering polar angle θ s is given in each panel.All distributions have been normalized to a maximum of one and the data has been averaged for both azimuthal incident angles of 0 • and 180 • in the MD simulations.

Fig. 6
Fig.6Experimental (black) and computed (red) angular distributions of the scattered H-atoms at different incidence polar angles θ s and kinetic energies.The maxima of the distributions have been scaled to one.The data has been averaged for both azimuthal incident angles of 0 • and 180 • in the simulations.

Fig. 7
Fig.7Impact positions of the H-atoms at the surface for the trajectories with incidence conditions θ i = 40 • , φ i =0 • and an initial kinetic energy of 1.92 eV.For reference the equilibrated surface is shown.The black lines highlight the surface unit cell, darker grey circles are oxygen atoms, lighter grey circles are Al atoms.The brighter the color of the points, the higher is the density of impact sites.

Fig. 8
Fig. 8 Model of the surface geometry resulting in different scattering angles for elastic scattering in case of an incidence polar angle of θ i = 40 • .

Table 1
Calculated and experimental lattice constants a and c and bulk modulus B 0 of bulk α-Al 2 O 3 87 .

Table 3
Sticking probabilities P sticking obtained in the MD simulations for different incidence conditions.

Table 4
Relative average energy losses obtained in MD simulations (sim.) and experiment (exp.)compared to the Baule limit (BL) for θ i = 40 • and θ s = 50 • .