An evaluation of the kinematic approximation in helium atom scattering using wavepacket calculations

We use 2-D wavepacket calculations to examine the scattering of helium atoms from dynamic assemblies of surface adsorbates, and in particular to explore the validity of the widely used kinematic scattering approximation. The wavepacket calculations give exact results for quasi-elastic scattering that are closely analogous to time-of-flight (TOF) experiments and they are analysed as such. A scattering potential is chosen to represent 8 meV helium atoms scattering from sodium atoms adsorbed on a Cu(001) surface and the adsorbates in the model move according to an independent Langevin equation. The energy broadening in the quasi-elastic scattering is obtained as a function of parallel momentum transfer and compared with the corresponding results using the kinematic scattering approximation. Under most circumstances the kinematic approximation and the more accurate wavepacket method are in good agreement; however, there are cases where the two methods give different results. We relate these differences to pathological features in the scattering form-factor.


Introduction
The technique of helium atom scattering has been recognised for many years as being particularly sensitive to the presence of surface adsorbates [1,2]. The early experiments of Karl Heinz Rieder and others [3][4][5][6] demonstrated that it is possible to determine adsorbate structure even for atoms, such as atomic hydrogen, that are difficult to observe using conventional diffraction techniques. The same experiments also hinted at the difficulties that would be involved in performing a fully quantitative analysis of helium scattering. In the early experiments, it was possible to avoid many of the difficulties using a hard-wall potential and a single scattering approximation for the dynamics. However, those approximations are only valid in the case of weakly corrugated systems.
For defects, such as surface-steps and isolated adsorbates, the scattering cannot be fully described without the inclusion of multiple scattering and the hard-wall model is difficult to justify at a quantitative level [7,8]. Thus, exact modelling of helium scattering from surfaces requires a fully quantum approach, with a realistic, soft, helium-surface potential. Such calculations are computationally expensive, even for fixed geometries, and only recently have they been used widely [8][9][10][11][12][13][14][15] Scattering from dynamical systems, such as in the study of surface diffusion using quasi-elastic helium scattering (QHAS) [16][17][18] represents a substantially greater challenge for calculation since the interaction potential itself is time varying. In addition, the results need to be averaged over an appropriate ensemble before results can be compared to experiment. In lieu of a rigorous solution to these problems, experimentalists have fallen back on the simplest approach, which is to use kinematic scattering [19]. In this approximation, each adsorbate is treated as a point scatterer as is the case of neutron scattering. It follows that many of the results derived in neutron scattering [20], and the related Fourier methods, can be carried over to the analysis of helium scattering. For example, with a moving adsorbate, or an assembly of moving adsorbates, the scattering amplitude is a function of time and the relevant ensemble averages, such as the intermediate-scatteringfunction, can be calculated simply using the convolution theorem [21,22].
Several recent measurements based on the Helium Spin-Echo (HeSE) method have led to a deep understanding of aspects of the adsorbate dynamics, such as: adsorbate-substrate interactions [23][24][25]; interadsorbate interactions [25][26][27]; adsorbate friction [22,28,29] and entropic contributions to surface mobility [30]. All of these works make use of an interpretive model that is based on a description of the adsorbate dynamics through a Langevin equation, coupled with an analysis of the scattering using the kinematic approximation. For better, or worse, the method has become the de-facto standard.
In the present work we analyse the validity of the kinematic scattering approximation in the study of adsorbate dynamics. Our approach is to use a fully quantum description of a helium wavepacket scattering from a dynamic assembly of adsorbates. We include a soft-wall potential, deduced from existing calculations of the surface chargedensity, and one that is able to represent the overlap of adsorbates. Thus all the main features of the experiment, such as the multiplescattering and the velocity spread in the beam are included. As before, we use a Langevin simulation to describe the adsorbate dynamics and the conventional kinematic approximation from the same adsorbate simulation gives a result for comparison with the fully quantum treatment. In order to make the calculations tractable and, in particular, to allow ensemble averages over times that are much longer than the mean hopping time, we have restricted the geometry of the 3-D dynamical system to 2-D. Thus, one dimension, z, is perpendicular the surface and the other dimension, x, is parallel to the surface. We have chosen parameters that model the helium scattering under typical experimental conditions, while the adsorbate motion is based on a 3-D model of Na adatoms moving on Cu(001). The calculations necessarily model a specific scattering-geometry and a particular system of interactions; however, we argue that the results have a greater generality and validity for the present methodology, based on kinematic scattering.

Methods
In order to represent quasi-elastic scattering using a realistic scattering method, we propagate a wavepacket of incoming helium atoms as it scatters from a time-dependent potential for the helium-surface interaction. The potential arises from the adsorbates as they move on the surface, and so its corrugated form fluctuates with time. As the wavepacket scatters from the fluctuating potential, it accumulates small energy changes. Averaging the outgoing intensity distribution over many wavepackets yields the quasi-elastic energy broadening, which we extract from the results in the same way that a time-of flight experiment is analysed [21,31]. The adsorbate motion is determined separately and is independent of the wavepacket dynamics. Thus, the subtle effects related to recoil of the adsorbate atoms during interaction with the helium [32] are not included in our analysis.
A complete simulation requires several components: (i) a description of the time-dependent potential that controls the scattering, given the instantaneous position of all the surface atoms; (ii) a prescription for the motion of the adsorbate atoms on the surface; and, (iii) a suitable wavepacket propagation routine. We elaborate on each of these points in the sections below.

Time-dependent scattering potential
The soft interaction and spatially extended nature of surface-interaction is one of the key characteristics that differentiates helium-from neutron-scattering. Our aim is to model those aspects of the repulsive interaction accurately. The repulsive part of the helium-surface interaction arises predominantly from the surface charge density and its dependence on the distance from the surface, z, can be approximated as an exponential function where V o and the stiffness, γ, can be extracted from calculations of electronic charge density [33]. z o is taken as a convenient reference plane. In the case of an isolated adsorbate there is a lateral corrugation of the charge density. Thus, the potential, as a function of the lateral coordinate, x, and the normal coordinate, z, becomes .
We find that, for an isolated adsorbate located at = x x , a the lateral corrugation is well approximated by a sinusoidal function where h is the corrugation amplitude and the width parameter, w, describes the lateral extent of the atom. For a moving adsorbate, so that the potential, as a whole, becomes time-dependent. Calculations of the surface charge in the appropriate density regime are available for Na atoms on a Cu(001) surface, having a range of lateral separations [34]. From these results we have extracted the parameters shown in Table 1 and, for an isolated adsorbate, the classical turning point of the resulting potential is shown in Fig. 1(a).
In the case of closely spaced adsorbates, overlap is handled by addition of the interaction for each adsorbate, in the region of overlap. Fig. 1(b) and (c) illustrates the form of the potential through the classical turning points for Na atoms at nearest-neighbour and next-nearestneighbour positions, respectively. We find that the contours in Fig. 1 Table 1 Parameters used to model the interaction of helium with Na adsorbates on a Cu surface. Note that the position of the reference plane z 0 is arbitrary, given a particular V t , and we choose a value that ensures the highest energies in the wavepacket are fully reflected at the surface.  agree with the trends demonstrated by Trioni et al. [34].
There is an important numerical advantage in limiting the lateral extent of an adsorbate corrugation to a width w. In the case of a moving adsorbate, the potential needs to be recalculated at each step of the dynamical simulation. To recalculate the potential in the whole scattering volume is a substantial task. By restricting the lateral extent of an adsorbate the volume that needs recalculation at each step is reduced significantly. For the same reason, we have not incorporated a long range attractive potential in our model. The result is a potential that describes the most important factors such as softness, adsorbate overlap and multiple scattering, within a computationally efficient prescription.

Adsorbate dynamics
The adsorbate dynamics are modelled using a Langevin equation where m i is the mass of the ith atom, R i is its position and V(R) is the "frozen" adsorbate-substrate potential. For the present work, the dynamical variables R i are one-dimensional, with motion taking place in the surface plane (hereafter notated as x). The exchange of energy between the adsorbate and substrate is determined by the friction coefficient, η, and the white-noise impulses, ξ i (t), are scaled according to the surface temperature, T, The term involving − F R R ( ) j i allows for a sum of pairwise interactions between the adsorbate atoms.
Following Alexandrowicz et al. [35], we use a sinusoidally modulated adsorbate-substrate potential and the same dipole-dipole repulsion model for the interaction between adsorbates, so that x Å is the distance between adsorption sites, = V 75 0 meV, = η 0.5 ps −1 and = p 0.58 Åe. The surface temperature is chosen to generate frequent hopping during the time taken for a single wavepacket to scatter from the surface. In this way the scattering dynamics samples the adsorbate dynamics efficiently. Trajectories that give the positions of all adsorbate atoms at every time-step and were generated by integrating Eq. (3) using a 4th order Runge-Kutta method, with a timestep of 0.5 fs and a surface temperature of = T 700 K. The trajectories are used in the calculation of the scattering-potential, Eqs. (1) and (2), as well for kinematic calculations, which we use for comparative purposes (see below).

Wavepacket propagation
For wavepacket propagation, we use the method of Tal-Ezer and Kosloff [36], along with the Fourier propagation algorithm of Kosloff and Kosloff [37]. The wavefunction, ψ(r, t), is represented on a discrete, 2-dimensional spatial domain which includes a region of free space, where = V 0, and a region containing the dynamic potential described above. The wavepacket starts entirely within the region of free space and is propagated in time using the operator,  U , where with and  H is the Hamiltonian operator, m is the mass of the atom described by the wavepacket and V(r, t) is the atom-surface potential. We ap- , such that the time evolution of ψ is then given by The quantity − ∇ ψ 2 is constructed by Fourier transforming ψ, multiplying by k 2 , where k is the corresponding wavevector, then inverse transforming. To reduce numerical errors, integration of Eq. (10) is carried out using the 4th order Runge-Kutta method. A typical calculation is illustrated in Fig. 2, which shows a series of snapshots of the wavefunction taken before, during and after scattering from the surface. The wavefunction is plotted logarithmically, as ψ*ψ, and is shown as a function of z in the left hand panels, and as a function of momentum, k z , in the right hand panels. Each panel, Fig. 2(a)-(e) corresponds to the same value of the lateral variable, x, and they cover a time difference of 29 ps, which is the time taken for the wavepacket to scatter completely. Unitarity was monitored throughout the calculation to ensure where = − ϵ 10 6 . If the unitarity condition was breached, the simulation was restarted with the timestep divided by two.
The amplitude of the wavefunction at the start of the calculation was defined by a window function, in z, chosen to minimise the range of k z present in the wavepacket. We use a Gaussian profile multiplied by a Hann window [38] so that where L, σ and z i are chosen to confine the initial wavepacket in a region away from the surface, where = V x z ( , ) 0 (see Fig. 2). For the simulations presented below we used a discrete grid of 512 × 512 pixels, corresponding to a (204.8 × 204.8)Å cell, with a typical propagation timestep of 1 fs and window parameters of = L 90 Å, with = σ 54 Å. Each wavepacket is propagated for approximately 29 ps, i.e. the centre of the wavepacket takes 14.5 ps to reach the surface and the same amount of time to return to approximately its starting position. The incident momentum of the wavepacket, k, was adjusted to represent 3 He with a mean beam energy of 8 meV and the window function (Eq. (12)) results in an energy spread of 0.1 meV. The energy spread is sufficient to resolve energy broadenings due to diffusive hopping of the adsorbate atoms. However, the length of the simulation time is insufficient to give an accurate ensemble average from a single calculation. The results below arise from a statistical average of 200 calculations. . 3 shows typical results generated for a single calculation described above. In the figure, the scattered intensity is presented on a grey-scale plot, as a function of the outgoing wavevector components, k ox and k oz . Fig. 3(a) shows the results for scattering from a single, static, adsorbate, using the potential shown in Fig. 1(a). A static adsorbate gives rise to scattering that is purely elastic as shown in the figure, where the main feature is the ring of scattered intensity lying on the Ewald sphere at = k o 3.4 Å − , 1 equal to the initial momentum of the wavepacket. Fig. 3(a) demonstrates that scattering occurs over a large range of final angles; however, there are significant variations of intensity with angle that arise from the scattering form-factor of the adsorbate. Fig. 3(b) shows results of scattering from a dynamical assembly of adsorbates. Specifically, 16 adsorbates, moving according to the Langevin model described above. The coverage, = θ 0.2, corresponds on average to 1 adsorbate every 5 adsorption sites. In this case, the motion of the adsorbates leads to energy transfers and the resulting intensity distribution is much more complex. The Ewald sphere, corresponding to elastic scattering is visible at = k 3.4 Å −1 and is accompanied by areas of inelastic scattering and a further ring at = k 4.8 Å −1 . The additional ring corresponds to inelastic energy transfer from vibrational motion parallel to the surface around the adsorption site. We also note that the elastic peak is wider than in the static simulation, corresponding to the expected quasi-elastic energy broadening.

Fig
The quasi-elastic broadening in data like Fig. 3(b) is extracted using the same process as is used in time-of-flight QHAS experiments [21,31].
Here the line-profile is taken to be the energy distribution in the incident wavepacket broadened by convolution with a Lorentzian energy profile. The latter is the profile expected for hopping motion [39]. The full-width at half maximum of the Lorentzian, ΔE(ΔK), is adjusted to give the best description of the data. Fig. 4 plots values of ΔE obtained in this way, as a function of ΔK, using open diamonds. For comparative purposes the results of a kinematic analysis of the same adsorbate dynamics are shown as a solid line. The kinematic analysis follows the method used by Jardine et al. [31], where the kinematic lineshape is first broadened by convolution with the energy resolution and then analysed identically to the data from the wavepacket calculation. The aim is to minimise systematic effects in comparison between the two scattering calculations that could arise from the particularity of the method used to extract information from the data.
We now discuss the features observed in Fig. 4 and compare the quantum calculations with the kinematic approximation. Looking first at Fig. 4(a), we show the results for an adsorbate having a height parameter h=1.61 Å, which best represents the case of a realistic Na adatom on a Cu surface (see above). The variation of ΔE with ΔK is complicated and, although the adsorbate motion is mainly single hops between adjacent sites the data cannot be described by a simple jumpmodel [40], which would have a sin 2 dependence dipping to zero at = K Δ 0 and = = K π a Δ 2 / 2.45 Å − , 1 with a maximum near = = K π a Δ / 1.23 Å −1 . Two factors contribute most to deviations from the sin 2 model. First, there are contributions from ballistic motion and dephasing of the lateral vibrations [31]. The effects are emphasised by the high temperature, = T 700 K, used to create the adsorbate motion and since these contributions increase at large ΔK, they obscure any dip in ΔE(ΔK) near ΔK ≈ 2.45 Å −1 . Second, the repulsive interactions, Fig. 3. Plot of scattered intensity for 8 meV 3 He, coded in a grey-scale, as a function of outgoing wavevector components for an 8 meV wavepacket incident on (a) a single static adsorbate, and (b) 16 mobile adsorbates following the dynamical model described in the text. k ox and k oz are the surface parallel and perpendicular components of the wavevector, respectively. In both cases, the simulation cell is 204.8 Å long and the spacing between adsorption sites is 2.56 Å. In (a) the main feature is the intensity lying on the Ewald sphere due to diffuse elastic scattering by the adsorbate. In (b) the elastic peak is also present, this time broadened by quasi-elastic energy changes, and accompanied by a further ring due to energy gain from vibrational motion, parallel to the surface, in the sinusoidal well, Eq. (5), as well as a broad inelastic background. coupled with the low dimensionality of the simulation leads to strongly correlated motion [22,35], which generates pronounced maxima and minima, such as the maximum seen at ΔK ≈ 0.25 Å −1 . The fact that both the wavepacket and kinematic calculations shown in Fig. 4(a) demonstrate these effects is confirmation that they arise directly from the adsorbate dynamics and are not artefacts introduced by the particulars of the analysis.
Turning now to the similarities and differences in Fig. 4(a) between the wavepacket results (open diamonds) and the kinematic results (solid line), we note that over most values of ΔK there is excellent agreement. The close similarity between the two methods is the most important conclusion from the present work, as it justifies the kinematic methodology that is widely used in the analysis of experiment [22][23][24][25][26][27][28][29][30].
Differences in Fig. 4(a) between the wavepacket and kinematic results are evident through the anomalies in two areas: (i) values of ΔK ≥ 2.4 Å −1 ; and (ii) values near ΔK ≈ 1.1 Å −1 . The differences at large values of ΔK can be attributed to a scattering angle near grazing emergence and the fact that, in this region, the scattered intensity in the wavepacket calculation is small so that the quasi-elastic peak cannot reliably be separated from the inelastic background. In order to investigate the origin of the difference between the two data-sets in Fig. 4(a) near ΔK ≈ 1.1 Å − , 1 we repeated the wavepacket calculations using a different scattering potential for the adsorbate atoms. Fig. 4(b) was calculated for an adsorbate atom with height parameter reduced to = h 0.683Å, which makes a significant change to the scattering form factor. Comparison of Fig. 4(a) and (b) shows that the anomaly near ΔK ≈ 1.1 Å −1 is absent in the case of small h (Fig. 4(b)). Thus, we can attribute the anomaly near ΔK ≈ 1.1 Å −1 in Fig. 4(a) to the form factor for scattering. Note that the other anomaly, at large ΔK, remains largely unchanged since the form factor for both values of height, h, result in little scattering at large angles. Fig. 5 provides a more detailed explanation of the anomaly near ΔK ≈ 1.1 Å −1 for adsorbates, having = h 1.61 Å. The top panel, Fig. 5(a), is the same as the previous figure and the vertical bar indiciates the area of concern. Fig. 5(c) shows the form factor for scattering from an isolated adsorbate, extracted from the variation of elastically scattered intensity with ΔK in the calculation shown in 3(a). The form factor varies strongly with ΔK as has been seen in other hardwall and wavepacket calculations [7][8][9]. In particular there is a deep minimum near ΔK ≈ 1.1 Å −1 . The effect of the minimum is seen in Fig. 5(b), which shows the difference in energy, β, of the maximum in the quasi-elastic peak from the mean energy of the wavepacket. A systematic deviation from zero indicates a degree of misfit between the observed lineshape and the ideal lineshape. Deviations are evident in both the range of both anomalous regions in ΔK. We can conclude that discrepancies between the wavepacket and kinematic analysis arise only in regions where the form factor behaves pathologically, for example close to the dip deep at ΔK ≈ 1.1 Å −1 .
These conclusions are supported further by results obtained at higher coverage. Fig. 4(c) shows the situation when the number of adsorbates is doubled, corresponding to a relatively high coverage, = θ 0.4. The adsorbate dynamics is now dominated by repulsive interactions and there are significant regions of overlap between adsorbates in the scattering potential. The resulting plot of ΔE(ΔK), Fig. 4(c), is correspondingly complex and as before, there is good agreement between the kinematic and wavepacket results. Interestingly the anomaly between the two methods, near ΔK ≈ 1.1 Å −1 and discussed above ( Fig. 4(a)), is much reduced. We attribute this to the increased overlap between the adsorbate potentials, and corresponding changes to the scattering form-factor, which will reduce the effect of the pathological     Fig. 3(a).
Surface Science xxx (xxxx) xxx-xxx dip noted in the discussion of Fig. 5(c). In fact, at these coverages no single form factor fully describes the scattering, yet we see that the linewidth derived from a kinematic analysis gives an accurate measure of the quasi-elastic broadening and one that is in good agreement with the full quantum treatment of scattering.

Summary & conclusions
Quasi-elastic scattering of helium atoms, and particularly the helium spin-echo method, has become an important approach for studying the dynamics of adsorbate motion on surfaces. In such work, the kinematic approximation is invariably used as the method to construct a model of the experiment, so that quantitative information on the true surface dynamics can be extracted by comparison with behaviour generated using Langevin simulations. In the present work we have presented the first quantum mechanical calculations of quasielastic scattering from a dynamic target. The scattering calculations are exact within the 2-D geometry we employ and a realistic soft wall potential has been constructed, based on recent ab-initio calculations for Na on Cu(100). Wherever possible, parameters for the various stages in the calculation have been chosen to best represent typical experimental conditions.
The approximations that typically lead to kinematic scattering [19] are the assumption of single-scattering from a weak, localised potential. In contrast, helium scattering from surfaces is best described as strong scattering from the delocalised potential of the surface electrons [1,41,42]. It is, therefore, surprising that the kinematic approximation and our exact quantum calculations are quantitatively equivalent over such a wide range of momentum transfer, ℏΔK. Not only are the methods comparable, the limited regions of difference between the two methods can be understood in terms of pathological features of the adsorbate form factor and, in particular, the strong dips that arises from multiple-scattering. It is unlikely that the agreement between results derived using kinematic scattering and those from exact calculations arises from a particular choice of interaction potential. Thus, it is fair to expect that our conclusions will have general applicability and, for example, the absence of an attractive part to the interaction potential is not in itself a significant restriction. The same argument can be applied to the dimensionality of scattering geometry, which would indicate that similar results should be expected in 3-D; however, such speculation will only be confirmed when calculations of sufficient complexity have been performed. Our conclusions provide strong support for the established methods of QHAS data analysis, at least in the case of atomic adsorbates.