Diffraction of fast heavy noble gas atoms, Ar, Kr and Xe on a LiF(001) surface, changing the tip of a ‘perfect’ AFM

Abstract

We investigate experimentally the diffraction of fast atoms of noble gas on a LiF(001) crystal oriented along the [100] and [110] directions. All exhibit some quantum features but wavelengths are so short that these effects are qualitatively described by semi-classical models. With increasing mass and energy, the scattering profiles show an increasing number of diffraction peaks forming an increasing number of supernumerary rainbow peaks but progressively weakening in contrast with the innermost peaks correspond to individual quasi-specular Bragg peaks disappearing first. Along the [100] direction, all observed azimuthal profiles are well described by Bessel functions allowing a simple semi-quantitative analysis. After removing the contributions of the attractive forces, we show how the surface corrugation amplitude and its variation strongly depend on the probing atom. These data should be compared with those accessible with an atomic force microscope (AFM).

---

The diffraction of fast atoms on crystal surfaces at grazing angle of incidence (GIFAD) has emerged in the last two decades as a new tool to investigate pristine crystal surfaces and to track the growth of thin films. It uses atoms with a kinetic energy E₀ ∼ keV and the same grazing geometry as reflection high energy electron diffraction, leaving the volume above the surface free for evaporation cells. In this grazing geometry, the perfect decoupling of the fast motion parallel to the probed crystal axis and the slower one perpendicular to it has been demonstrated for elastic diffraction.^1–3 By tuning the angle of incidence θ_i between typically 0.2° and 2°, the effective interaction energy E_⊥ = E₀sin²θ_i can be adjusted from a few meV to several eV. In the elastic regime, GIFAD is equivalent to thermal energy helium scattering (TEAS) but it differs in the inelastic regime because the momentum transfer to the surface atoms occurs in several gentle collisions so that elastic diffraction can be observed at larger values of E_⊥, higher surface temperature and with heavier projectiles. We investigate here experimentally and via simulations the consequences of increasing the mass of the projectile atoms.

• For comparable energy the wavelength reduces (see Table 1) allowing potentially higher spatial resolution but the Bragg angle ϕ_B reduces with λ requiring more collimated beams.

Table 1 De Broglie wavelength λ in pm for different energies and atoms. On LiF(001), the distance a_⊥ between identical atomic rows along [100] and [110] are 202 and 286 pm respectively corresponding to reciprocal lattice vectors G_⊥ = 3.11 and 2.20 Å⁻¹. The associated Bragg angle is ϕ_B ≃ λ/a_⊥

Energy	He	Ne	Ar	Kr	Xe
1 keV	0.45	0.20	0.14	0.098	0.079
100 meV	45	20	14	9.8	7.9
10 meV	144	64	45	31	25

---

• The momentum transferred to the surface increases rapidly, increasing significantly the probability of phonon excitation as described by the Debye–Waller factor adapted to GIFAD.^4–6

• The larger number of valence electrons increases the magnitude of the Pauli repulsion, pushing the classical turning point away from the surface. This is balanced by the increased polarizability attracting the projectile towards the surface.

Within the Born–Oppenheimer approximation, the diffracted intensities correspond to the quantum scattering of the projectile atom in the potential energy landscape (PEL) describing the energy of the projectile at any location E(x,y,z) above a surface lattice unit. A complete approach such as e.g. wavepacket propagation^7–9 or close-coupling¹⁰ provides exact correspondence between the PEL and the diffraction patterns at various energies, as well as possible bound state resonances^11,12 or quantum reflection.¹³ This is specific to the quantum regime where λ is larger than typical dimensions of the surface topology. The first direct consequence of the reduced wavelength attached to larger masses is that the observations presented here correspond to the semi-classical regime and can, at least qualitatively, be explained by phases attached to classical trajectories and interferences between different paths ending at identical final scattering angles θ_f, ϕ_f. In this regime the connection to the surface topology is more direct and the information derived should compare with that from atomic force microscopy (AFM). For example, the corrugation amplitude describing the size of bumps measured by an AFM on top of the surface is also measured by GIFAD but in the reciprocal space. Such analogy was popularized by the semi-quantitative hard corrugated wall model (HCW).¹⁴ It is an optical model considering straight line reflection from the iso-energy surface z_E(x,y) assumed to represent the surface where the classical trajectories are reflected. The diffracted intensities are therefore obtained by a Fourier transform of the topology. However, as a reciprocal space technique, atomic diffraction measures scattering angles acquired essentially while bouncing off these bumps. However, trajectories are also affected by angular changes due to the attractive forces modifying outgoing trajectories.

The vast majority of publications on atomic diffraction on surfaces consider He as a projectile (see e.g. ref. 15 for a recent review at thermal energies) but numerous theoretical studies^16–19 as well as a few experimental works considering Ne.^20–24 Investigations using Ar on different surfaces are scarce with early attempts to measure physisorption and trapping via mass balance techniques²⁵ then, in the eV regime^26–28 scattering profiles somewhat similar to Fig. 1 were recorded point by point revealing classical rainbows and, together with energy loss measurements indicated early signs of transitions from binary collisions to the axial channeling regime²⁹ which is fully established in GIFAD.^1,2,7,30 With Kr and Xe it seems that only classical descriptions of the projectile were involved.^31–33

Fig. 1 Raw camera capture of the scattered atoms on the detector for 500 eV Ar impinging LiF(001) along [110] direction at θ_i = 1°, ϕ_i = 0°. The dashed circle θ_f = θ_i corresponds to energy conservation; |k_f| = |k_i|, the vertical line is the specular plane k_yf = k_yi. In the lab frame (inset), the azimuthal or lateral deflection angle is ϕ_f = arctan(k_yf/k_xf) ≃ k_yf/|k_i|, less than a degree. In the (y,z) detector plane, the primary beam arrives from top with k_zi = k₀ sin θ_i ≡ k_p and the angle Θ_f = arctan(k_yf/k_zf) ≃ arctan(ϕ_f/θ_i) is equivalent to that measured in TEAS and can reach ±90°.

Detailed data for helium^30,34 and neon^24,35,36 have already been published over a wide range of energy and angle and only the new results with Ar, Kr and Xe will be presented. Rather than three independent sections with each projectile, the paper focuses on specific topics and their evolution along the Ar, Kr and Xe sequence. The outline is as follows: after a brief recall of the GIFAD technique in Sections 1, Sections 2–6 are a purely experimental description of the evolution from quantum to classical features. The elastic and inelastic diffraction peaks are presented in Sections 2 and 3 respectively. Section 4 presents the supernumerary rainbows while Section 5 focuses on the classical rainbow and Section 6 focuses on comparatively larger collision energy where no trace of quantum behavior can be seen but where specific features can be efficiently described by semi-classical models.

Section 7 presents all the data in a more comprehensive way highlighting the role of the attractive forces extending up to eV effective energies and that of the topology, dominant at larger energies. A model potential energy landscape fitted to the diffraction data is presented for the [100] direction and the resulting surface topology is derived for all noble gases along this direction before tracing perspectives and conclusion in Section 8.

1 Experimental setup

A detailed description of the GIFAD setup and of the associated coordinate system can be found in ref. 37 and 38 respectively and only a simple sketch is presented below. An ion beam of kinetic energy E₀ in the keV range is neutralized in a charge exchange cell and severely collimated by two diaphragms half a meter apart to reduce the angular divergence to the 0.1 mrad range. It impinges on the crystal surface at a grazing angle of incidence θ_i ∼  1°, usually along a low index crystallographic direction, and all diffracted beams are within a narrow cone with an opening angle close to θ_i. A microchannel plate imaging detector³⁹ located ∼1 m downstream records all diffracted beams at once. When elastic diffraction is present, sharp spots are visible on the energy conservation circle (|k_f| = |k_i|) while inelastic diffraction occurs more or less close to this circle depending on the surface temperature.^3,40

Fig. 1 displays a scattering pattern typical of inelastic diffraction by heavy projectile with short wavelength where mainly supernumerary rainbow structures are visible as discussed in Section 4. No sharp spot is present and therefore the specular circle is not directly visible. It also illustrates the angular coordinates used here, θ and ϕ are referred to as the polar and azimuthal angles in the lab frame. The azimuthal (lateral) deflection will be reported in degree or in multiple values of the Bragg angle. When comparing different energies or different incidence angles we use the relative angle Θ_f = arctan k_yf/k_zf ≡ arctan ϕ_f/θ_i to suppress the trivial dependence of the azimuthal scattering ϕ_f with the angle of incidence θ_i. This angle then compares directly with the coordinates of a TEAS setup where the incidence angle is referred to the surface normal. Assuming that the parallel velocity can be ignored, Fig. 1 would represent the scattering of a 150 meV Ar beam directed exactly normal to the surface. However, the interaction potential would be that of the actual surface but averaged along the parallel velocity. Each atom in the lattice cell is replaced by a translation invariant row of aligned atoms. Such rows are separated by a distance a_⊥ producing the Bragg angle ϕ_B and their linear density is 1/a_‖ such that the surface of the lattice cell a_⊥ × a_‖ is constant. The coordinate Θ is also more appropriate to describe the refraction effect relating the deflection angle Θ_surf acquired when bouncing along the sides of the atomic rows to the one Θ_obs modified by the attractive forces when escaping. This later is brought closer to the surface and the relation is well-described by the Snell–Descarte law:

where D is the well-depth resulting from attractive forces so that E_⊥ + D represents the effective energy when hitting the surface and the velocity or refraction index entering Snell's formula.

2 Tracking the elastic diffraction

Elastic diffraction is identified as a narrow peak at the specular angle in the polar scattering profile (see e.g. Fig. 3 in ref. 6). Each elastic diffraction spot has a profile very close to that of the primary beam, the latter being limited by the size of diaphragms and detector resolution. For He and Ne projectiles on the LiF(001) surface, the ratio of elastic to inelastic diffraction was found independent of the crystal direction and could be modeled by the Debye–Waller factor adapted to GIFAD.^4,5 At room temperature the following scaling was observed⁴⁰ with M_p the projectile mass in Dalton, E₀ its total kinetic energy in meV and θ_i the incidence angle. For E₀ = 1 keV and θ = 0.5°, the exponential term amounts to 70, 19, 3.6, 0.08 and 0.002% for He, Ne, Ar, Kr and Xe respectively. On top of this formula, an empirical pre-factor A(M_p) of 0.55 for He and 0.16 for Ne was observed⁶ probably due to our limited surface coherence and indicating that the sensitivity to these defects increases with the projectile mass. When analyzing the polar profile reported in the inset of Fig. 2, we measured a very weak elastic contribution representing only 0.3 ± 0.3% of the scattered intensity. Note that both the Gaussian peak and the log-normal peaks are sitting at the specular scattering angle. They appear shifted because the median value of a log-normal is not at maximum. The narrow Gaussian peak has the same width as the primary beam while the inelastic contribution is almost ten times higher and twenty times broader. This weak elastic ratio can also be identified but less quantitatively by using a double differential filter to the full scattering image, as detailed in ref. 41. This elastic ratio indicates a pre-factor below 10% for the effective Debye–Waller factor with Ar on our LiF crystal, consistent with the one of 55% and 18% measured with helium and neon projectiles on similar samples as detailed in Fig. 6 of ref. 6. This suggests that larger elastic diffraction could probably be observed at lower surface temperature.⁶ Our samples were cleaved by hand after irradiation (see acknowledgements) and a more elaborate procedure could also improve the surface quality. All data presented here were recorded at room temperature and, for each data set, the surface quality was checked by helium diffraction. No sign of elastic diffraction is identified for Kr and Xe.

Fig. 2 For 500 eV Ar impinging LiF along [100] at θ_i = 0.29°, the diffracted intensity measured on the specular circle is fitted (red line) by quasi Lorentz line-shape (eqn (2)) measured along the [Rnd] direction. The green lines are for the individual contribution and the measured intensities are reported on the left. The inset shows the polar profile with a weak Gaussian peak at specular angle θ_s = 2θ_i on top of a broad log-normal profile⁶ at the same location. The arrow points to the expected intercept of the sample surface plane at θ_i.

3 Tracking inelastic diffraction peaks

Along the azimuthal direction, inelastic peaks can be well identified as long as their line-width σ_lw is less than the peak separation i.e. the Bragg angle ϕ_B.

3.1 The inelastic azimuthal line profile

As proposed in ref. 42 and 43, the inelastic azimuthal line profile does not depend of the crystal orientation and can easily be measured by orienting the surface far away from a low index direction also called random direction [Rnd] where only one diffraction peak is present. This was investigated in detail in the quasi elastic regime with He and Ne projectiles where a quasi-Lorentzian line shape with reduced wings was proposed⁴³ in the form of the product of a Lorentzian by a Gaussian

LG(ϕ) = A·e^−ϕ2/2w2/(4w² + ϕ²), (2)

where the standard deviation σ_lw is related to the width parameter w by σ_lw ≃ 0.7w and A is a normalization factor. All the line shapes used here and plotted in green in Fig. 2–5, use this form with a width parameter interpolated from measurement along the [Rnd] direction⁴⁴ and referred to by its standard deviation σ_lw.

Fig. 3 Along the [100] direction, scattering profiles of (a) 4 keV Kr and (b) Xe atom at θ_i = 0.45° and 0.48° respectively, have been fitted by a HCW formula forcing the relative intensities to follow I_m = J_m²(ζ). Both profiles give comparable values of ζ close to 3.9 for which the intensity I_±1 = J₁²(ζ) is close to zero (see text). The linewidth σ_lw exceeds the Bragg angle ϕ_B. Data for Xe have been symmetrized.

Fig. 4 Scattering profile of 700 eV Ar impinging LiF along [110] at θ_i = 0.42° (E_⊥ = 38 meV). It is fitted (red line) by line-shape measured along a [Rnd] direction and no constraint on diffracted intensities I_m. The peaks labeled s_i correspond essentially to even diffraction orders |m| = 0, 2, 4 while |m| = 7, 8, 9 contribute to the outer classical rainbow.

Fig. 5 Scattering profile of 1784 eV Ar along [110] at θ_i = 0.5°. The red line is a fit with a line-shape width σ_lw ≈ 30% larger than the Bragg angle and no constraint on the intensities. Close to the center, the alternance bright/dark of even/odd order persist and supernumerary peaks are dominated by a single diffraction order while more and more contribute to the outer peaks. The insert shows a lower incidence of 0.43° with less supernumerary peaks but more contrasted ad almost identical to the value displayed in Fig. 6 at 500 eV and θ_i = 0.83° and E_⊥ = 104 meV.

3.2 The [100] direction

Fig. 2 shows an azimuthal scattering pattern corresponding to a narrow band around the specular circle⁴⁵ for 500 eV Ar atoms. The inelastic peaks are resolved but overlap significantly. Using the line profile recorded along the [Rnd] direction, the only free parameters of the fit are the diffracted intensities I_m reported in the left-hand side. These intensities are very close to the values I_m = J_m²(ζ = 1.52) where J_m are the Bessel functions of order m so that the full diffraction pattern could have been fitted by a single number ζ = 1.52. This nice feature corresponds to the predictions of semi-classical perturbation theory² as well as that of the hard corrugated wall model (HCW). In the latter, the PEL above the surface along [100] is represented by equipotential lines with a pure sinusoidal shape providing a simple interpretation for ζ = k_⊥z_c and z_c is called the corrugation amplitude. For instance the value of ζ = 1.52 derived above translates into z_c ∼ 10 pm for the associated wave vector k_⊥ = 15.6 Å⁻¹. We will see in Section 7 that this value is almost twice too large while the interpretation in terms of rainbow angle seems more correct. Using geometric optics, the deflection angle derives from the inclination (blaze) of the corrugation function Θ = −2 arctan ż(y) and the classical rainbow angle Θ_rb is defined as the maximum deflection corresponding here to y = ±a_⊥/4:

The formula offers a simple estimate of the rainbow angle in conditions where only a few diffraction orders are present as in Fig. 2 for which the rainbow is not visible as a distinct feature. Along the [100] direction, the same agreement with Bessel functions was observed for He and Ne projectiles (see Fig. 3 of ref. 7, Fig. 3.12 of ref. 30 and Fig. 9 of ref. 36).

For Kr and Xe, we measured, along the [Rnd] direction, an inelastic width σ_lw larger than the Bragg angle and, consistently, we could not clearly resolve adjacent diffraction peaks. However, one can easily identify situations where the m = ±1 diffraction peaks are weak because their intensity oscillates in quadrature with m = 0 and m = ±2, offering the opportunity to identify diffraction peaks in spite of σ_lw > ϕ_B. The first zero in the J₁(ζ) Bessel function occurs for ζ = π + π/4 ∼ 3.9 where the π/4 is interpreted as the Gouy phase or Masslov index ‘naturally’ present in the Bessel function. The scattering profiles displayed in Fig. 3 are close to this criterion.

3.3 The [110] direction

Fig. 4 shows a scattering profile recorded of Ar atoms along the [110] direction, together with a free fit of the intensities using a line-profile fitted along the [Rnd] direction. Once again, peaks are identified mainly because the intensity of the even and odd diffraction orders oscillate in quadrature in the paraxial region of small deflection³ (see e.g. Fig. 2c and 3c of ref. 6 for illustration).

Along this direction, as already noted with He and Ne,^24,34 data modeling with a HCW requires a second term in the Fourier expansion of z_E(y) to describe the additional electron density due to Li rows in between the F rows. The additional parameters make the fit less stable and the analysis more complex as discussed at the end of Section 7.1). The HCW was not used along [110].

4 Tracking the supernumerary rainbows

Another direct evidence of quantum behavior is the presence of supernumerary rainbows. Their observation in natural conditions was essential in the foundation of the wave nature of light⁴⁶ and are quantitatively described by the Mie theory of light scattering.⁴⁷ More generally, supernumerary rainbows are ubiquitous in quantum scattering theory. A classical rainbow is associated with an extrema in the deflection function and, most often, on each side of this extrema, two different semi-classical trajectories are scattered at the same final angle and may interfere. The associated phase shift between these trajectories increases as these trajectories separate and supernumerary rainbows are observed when the phase shift reaches a multiple of 2π. Therefore, supernumerary rainbows indicate a situation with λ ≪ z_c where pure quantum behavior, such as bound state resonance^12,48 or quantum reflection^13,49 can probably be neglected. In GIFAD, increasing the number of diffraction orders rapidly leads to a quasi-continuous profile,^50,51 as if the discrete values of J_m²(ζ) would converge to the envelop J_ν²(ζ) where J_ν is the uniform Bessel function such that J_ν(ζ) = J_m(ζ) where ν = ϕ_f/ϕ_B (see e.g. Fig. 2.5 of ref. 30). This would be exact if the line profile were negligible, which partly contradicts the fact that individual peaks cannot be resolved so that this ideal situation is never reached.

4.1 Along the [110] direction

Along the [110] direction, such scattering profiles where individual diffraction peaks seem absent have been observed with helium,^34,51 neon^24,35 and argon (Fig. 1). Two intense classical rainbows are present at the maximum lateral deflection. In between, seven weaker maxima associated are present with four supernumerary rainbows for positive or negative deflection.

Fig. 5 shows a similar profile in a more quantitative plot where the supernumerary rainbows are numbered s₁ to s₄ starting from the outer classical rainbow position labeled as ±r_b. Fig. 5 also displays a free fit by individual lines showing that the inner supernumerary rainbow peaks s₂ to s₄ still correspond to individual diffraction peaks m = 0, ±2 and ±4 separated by 2ϕ_B as in Fig. 4.

This outlines the smooth continuity between scattering patterns with resolved inelastic peaks and supernumerary rainbows but also the associated difficulty because the position of the inner peaks is first attached to the location of actual diffraction peaks while it evolves continuously in the continuous Bessel function.

Fig. 6 reports the evolution of the scattering pattern P(Θ_f) during a θ-scan where the incidence angle θ_i is varied. Note that to avoid overlap between profiles, low incidence angles are plotted on top. The number of supernumerary rainbows increases from two to four between θ_i = 0.4° and 1.04°. However, the attenuation of the contrast makes their identification more hazardous while the classical rainbow peak becomes rapidly dominant.

Fig. 6 Θ scattering profiles recorded on the specular circle for 500 eV Ar along LiF[110]. The angle of incidence θ_i of the nine profiles is indicated on the left hand side and the associated value of E_⊥ is on the rhs.

The insert in Fig. 5 shows that scattering profiles recorded at different energies and angles but corresponding to identical values of E_⊥ are almost indistinguishable when plotted on a common Θ_f scale. The strict equivalence was established for the elastic intensity^2,10 but seems to apply for inelastic intensities on the specular circle θ_f = θ_i suggesting that the inelastic intensities are probably also close to the elastic one at θ_f = θ_i.

Similarly to ref. 35 for Ne, the angular location of the rainbow and supernumerary rainbow peaks is reported using the angle Θ = arctan(ϕ_f/θ_i), in Fig. 7 for Ar.

Fig. 7 Evolution of the rainbows peaks with E_⊥ for 1784 eV Ar along [110] direction. The positions are measured by hand and the continuous black line is a Snell transform of the horizontal dashed line with D = 50 meV.

For heavier projectiles at comparable energies, more supernumerary rainbows are present, but since the linewidth increases while the peak spacing decreases, the contrast rapidly weakens. An effect comparable to the contrast weakening visible in the bottom profiles in Fig. 6. For Kr, only three supernumerary rainbows could be unambiguously identified and only at low energy. These are reported in Fig. 8. For Xe, Fig. 9 shows that supernumerary rainbows might be present, but their identification is limited both by resolution and statistics.

Fig. 8 For 4 keV Kr atoms along the [110] direction, the peak position of the classical outer rainbow and inner supernumerary rainbows is plotted as a function of E_⊥. The full black line corresponds to an optical model using Snell's law with Θ₀ = 19.2° and a well-depth D = 92 meV. The inset shows the evolution of the rainbow angle for 1 keV Xe and the full line is for Snell's formula with Θ₀ = 13.2° and D = 135 meV.

Fig. 9 Azimuthal scattering profile of 1 keV Xe impinging LiF[110] at θ_i = 0.4° and 0.8°. The statistics is poor but indications of two supernumerary rainbows are present at θ_i =0.4° and possibly three or four at θ_i = 0.8°. The full line is a B-spline through the data points.

4.2 Along the [100] direction

Along the [100] direction and at low value of E_⊥, the corrugation amplitude is typically 3–4 times smaller than along the [110] and the number of supernumerary rainbows is reduced by the same amount. Disregarding the deconvolution in terms of diffraction orders, Fig. 3 shows the presence of a central maximum which is a direct evidence of a single, weakly contrasted supernumerary rainbow with Kr and Xe. For Ne, Fig. 8b of ref. 36 shows a better resolved but comparable profile while Fig. 10b indicates that, for Ar, a similar situation with I_±1 ≃ 0 (ζ ∼ 5π/4) is reached with E_⊥ ∼ 250 meV. More supernumerary rainbows probably develop at larger energy E_⊥ but, at room temperature σ_lw washes out the central structure where peaks are separated by 2ϕ_B.

Fig. 10 (a) Diffracted intensities I_±m = (I_−m + I_m)/2 for 500 eV Ar atoms along the [100] direction. Below 50 meV all intensities I_m were left free but above, a fit forcing a Bessel form I_m = J_m²(ζ) was used. (b) The derived rainbow angle () on the lhs is a good description of the actual value whereas, interpreting ζ as a corrugation amplitude z_c = ζ/k_⊥ on the rhs leads to an unphysical increase at low energy. Considering refraction with eqn (1) and D = 50 meV provides a more plausible value of the actual corrugation ().

5 Tracking the rainbow angle

Before the discovery of GIFAD, the rainbow angle was a reference measurement to compare theoretical PEL with scattering experiments under axial channeling conditions^18,52,53 with interaction energies E_⊥ above a few eV.²⁸ At these energies, the rainbow structure is the most salient feature of the scattering profile but tracking its position accurately is not straightforward. The rainbow has a classical origin but the exact shape is deeply affected by quantum features. In the quantum regime where λ > z_c as in Fig. 2 it does not correspond to any visible feature. Then, when supernumerary rainbows start to emerge, the shape of the rainbow peak becomes an Airy-like profile (see e.g. eqn (13) of ref. 30) affected by the inelastic line-width σ_lw associated to the thermal agitation of the surface atoms. At larger energy E_⊥ and larger values of σ_lw the rainbow profile is dominated by σ_lw and the point of maximum intensity is considered a fair estimate of the rainbow position. A simple alternative is to use the statistical width of the entire scattering profile, either from quantum point of view: σ_m² = Σ_mm²I_m⁵⁴ or from the classical one σ_ϕf² = Σ_i(αi)²I_i where I_i is the intensity measured on the channel i of the histogram of scattering angles normalized to unity and α the angular calibration. Both values should be close as they are connected: σ_ϕf² ∼ ϕ_B²σ_m² + σ_lw². This robust and parameter-free definition allows a smooth continuity through the different regimes as shown in Fig. 13 and Fig. 21 of ref. 34.

We use here the coordinate Θ and the angle Θ_rb which offers the double advantage to correspond to the deflection angle measured in TEAS and that to suppress the linear dependence with θ_i or k_iz (see e.g. insert in Fig. 5 or compare Fig. 11b) and Fig. 12 using Θ and ϕ scales respectively).

Fig. 11 Rainbow scattering angle of (a) Kr and (b) Xe atoms along the [100] direction of LiF[100]. The rainbow angle is tracked by fitting the scattering profile with Bessel functions as in Fig. 3. The straight dashed lines have been adjusted by hand so that the red line corrected by refraction in eqn (1), pass through the data. The well-depth used are D = 92 and D = 135 meV for Kr and Xe. Note that E_⊥ goes up to 16 eV.

Fig. 12 Raw 2D (θ_f,ϕ_f) scattering profiles of 4 keV Xe atoms (λ_‖ ∼ 40 fm) along [Rnd], [110], and [100] where ϕ_B is 0.008° and 0.011°, respectively. The square indicates the angular scale with the vertical arrow on its side referred to the surface plane θ_f = 0°. The energies E_⊥ (eV) are: 0.19, 1.84, 4.08, 11.2, 21.8 for [110], 0.14, 0.41, 1.24, 2.89, 6.67, 15.2, 22.9 for [Rnd] and 0.19, 0.99, 2.44, 4.57, 9.15, 16.2, 22.9 for [100].

5.1 Along the [110] direction

For all energies investigated here, the supernumerary rainbow is well developed along the [110] direction indicating λ ≫ z_c (or ζ ≫ 1) and the rainbow position was estimated by hand by trying to point the maximum of the outer peaks. The corresponding values are plotted in Fig. 7 and 8 for Ar and Kr, respectively. They show a stable value at large values of E_⊥ and a rapid increase at low energy. In these figures, the solid black line running through the data corresponds to a refraction model (eqn (1)) with a fixed emission angle Θ_surf of 19° and 25° and a well-depth D of 50 meV and 92 meV for Ar and Kr respectively. For 1 keV Xe with E_⊥ between 100 and 800 meV, the rainbow peak position is reported in the inset of Fig. 8 together with a form: with D = 135 meV.

5.2 Along the [100] direction

The rainbow angle is more difficult to point directly along the [100] direction, as it is typically Θ_rb ∼ 6°, three to four times smaller than along [110]. Taking advantage of the good fit with Bessel functions J_m²(ζ) in Fig. 2 and 3, we use this very stable data reduction procedure to extract the number which, within the HCW is the rainbow angle as long as Θ ∼ tan Θ. Note that with the same cosine shape and small-angle approximations Θ_rb is also directly connected to the standard deviation discussed at the beginning of the section, .

Fig. 10b and 11a, b show the measured rainbow angle for Ar, Kr and Xe respectively. The values are much smaller than those observed along the [110] direction in Fig. 6–8. As for the [110] orientation, the classical rainbow angles for the [100] orientation increases significantly for energies below 1 eV for all projectiles and this increase is well-modeled by a refraction model using the values listed in Table 2 close to theoretical values⁵⁵ of the Well-depth.

Table 2 Parameters used to model our data by the refraction effect in eqn (1): with α in . The energy range is limited along [110] and α_[110] was taken as 0 while the broader range investigated along [100] forces the use of an energy dependence. The corresponding corrugation amplitude is not z_c = G_⊥Θ_obs but z_c = G_⊥Θ_surf. D₀ are the theoretical values in ref. 55

	D 0 (meV)	D (meV)	Θ 0[100]	α [100]	Θ 0[110]
Ar	−62.5	−50	6.24°	5.2	23°
Kr	−94.2	−92	4.7°	4.1	19°
Xe	−152	−135	2.76°	3.2	13°

---

6 Larger impact energy, classical profile

We use the same HCW model in spite of the fact that all the trajectories should add up incoherently (without attached phase). This is not a problem since using the experimental line width of a few Bragg angles washes out all interferences converging to the classical envelope. In addition, the optical treatment of the slowly varying phase properly deals with the complex problem of rainbow divergence in classical scattering.

Fig. 12 displays the evolution of the raw 2D scattering profiles recorded for 4 keV Xe atoms along the [110], [Rnd] and [100] directions between a few tens of meV and a few tens of eV. Apart from the lowest energies hardly visible in the bottom but already plotted in Fig. 9 and Fig. 3b), most of the scattering profiles have no or very weak residual quantum structures but they still have distinct features such as the central spot on the top left side. The latter is characteristic of the [110] direction and is observed for all noble gases appearing at increasing values of E_⊥ along the He-Xe sequence. For He it is visible above 2 eV, see Fig. 2 and 3 of ref. 30 and Fig. 18 of ref. 34, while it appears here for Xe for E_⊥ 10 eV. It corresponds to a second rainbow peak which is neither a secondary nor a supernumerary rainbow and will be discussed in the next section. Fig. 12 shows that at large values of E_⊥ the rainbow angle associated with the [100] direction is significantly larger than that of the [110] direction. This is in contrast with the low energy range, E_⊥ < 0.5 eV, where the rainbow angle along [110] is close to 30° while it is closer to 5–8° along [100] for all noble gases. Note also that the sharp increase of the relative width σ_Θ ∼ σ_ϕ/θ_i visible in Fig. 11b) is hardly seen in the absolute scale of Fig. 12 where σ_ϕ is visible. Fig. 12 also shows that the line width defined by the scattering profile recorded along the [Rnd] direction becomes much larger than the Bragg angle ϕ_B and finally compares with the width of the rainbow peak.

7 Data analysis and modeling

Before discussing the minor differences along the He–Xe sequence, we analyze the common feature appearing in the high-energy range where the atomic arrangement along [100] or [110] are decisive and, in the low-energy range where attractive force plays the dominant role. Since no specific quantum effect has been observed, we use the simple HCW model already mentioned in the data reduction procedure.

The hard corrugated wall model considers free propagation of the atom and instant reflection on a corrugated wall. In terms of physics it acknowledges that, due to the exponential character of the Pauli repulsion, most of the momentum transfer is indeed located close to the equipotential surface of energy E_⊥. As an optical model it takes into account interferences between all rays (trajectories). This does not mean that constant velocity is a good approximation for the dynamics of the atoms, it simply means that above the surface, most of the phase differences arise from the path difference due to the topology described by the corrugation function. This was illustrated recently² using first order perturbation theory with a realistic Morse potential along z and shifted along y according to a cosine shape. The model has the full dynamics of the incoming trajectories inside the Morse potential but the diffracted intensities are still given by J_m²(k_⊥z_c). This is because the smooth acceleration and more sudden deceleration when approaching the iso-energy curve is the same for all incoming trajectories leaving the path difference due to the topology z(y) as the only contribution. The correction due to the outgoing trajectories appears with second order perturbation theory⁵⁶ and, placing the potential well at different location⁵⁷ also breaks the direct correspondence with the HCW.⁵⁸

7.1 High energy part, detailed topology

The higher the energy, the closer we approach the surface and the larger the contrast of the equipotential line because the region between surface atoms offers less resistance. This describes well the [100] direction where no surface atom lies between the atomic rows.

For He and Ne, it was already noted that the evolution of the corrugation amplitude is quite different along the [100] and [110] directions.^36,59 The simple interpretation is that along the [100] direction the surface appears as made of a single type of atomic row (see inset in Fig. 13). The maximum of the PEL sits on top of such a row where F⁻ and Li⁺ ions alternate while the minimum lies in between. At higher energy, the minimum decreases faster (is softer) than the maximum resulting in a corrugation amplitude that increases with energy as visible above E_⊥ = 0.5 eV on the rhs scale of Fig. 11 where the value ζ derived from the HCW is converted into a corrugation amplitude. Along the [110] direction the lowest point in the PEL corresponds to a row of Li⁺ ions, these have a reduced electron density but their K shell electrons are more deeply bound and more difficult to penetrate while the larger electron density of the F⁻ row is comparatively softer so that the corrugation amplitude is more or less constant at low energy (see e.g. Fig. 3–10 and 3-3 of ref. 30 resp.) with a tendency to decrease at energies E_⊥ above a few eV and in Fig. 12 for Xe.

Fig. 13 Real space corrugation z_c = z_E(y_top) − z_E(y_bot) of the average [100] potential energy landscape as a function of E_⊥ plotted in a scale. z_E(y) are the equipotential lines of the PEL fitted to the diffraction data using the Hardwall model with angular refraction (Beeby correction). The dotted lines correspond to energy domains where no experimental data were available. Data for He and Ne are from ref. 34 and 36 respectively.

This can be analyzed in more detail by investigating the exact shape of the PEL but the most convincing evidence of the presence of the Li⁺ row is the bright spot appearing on the top left of Fig. 12. As mentioned above the same spot has been observed for He, Ne, Ar and Kr at increasing values of E_⊥ and it has a simple classical origin in terms of a second rainbow structure. Assuming that the energy E_⊥ is large enough to reveal directly the presence of the Li⁺ row in the center of the lattice cell, i.e. that the iso-energy curve shows two maxima, a big one on the F⁻ row and a smaller one on the Li⁺ row, the deflection function will have two extrema associated with the two inflection points in the PEL and two new rainbow peaks at positive and negative values of the associated scattering angle. Before reaching this situation, the bottom of the PEL flattens and a peak appears in the scattering distribution and its intensity increases progressively.

For He at E_⊥ = 4.2 eV, this was modeled in Fig. 18 of ref. 34 using the HCW model with a corrugation function having two Fourier components Z(y) = h₁ cos ỹ + h₂ cos 2ỹ with ỹ = G_⊥y and h₂ ∼ h₁/10. For the pattern in Fig. 12 with Xe at E_⊥ ∼ 23 eV, the HCW using h₁ = 18.8 pm and h₂ = 3.6 pm produces a good fit made of ∼250 diffraction orders which are convoluted by a broad Gaussian profile with σ_lw ≃ 30ϕ_B. The width σ_lw is so large that the same result can certainly be derived from classical models without any interference but the HCW formula (eqn (2) in ref. 34) is so simple that the calculation is instantaneous and covers the full semi-classical range merging smoothly to classical behavior.

In practice though, for the values of h₁ and h₂ used the iso-energy curve Z(y) does not yet have a double maxima shape, however the bottom of the corrugation function Z(y) becomes so flat that the quasi specular reflection in this region produces the intense maximum observed in the center of the scattering profile. Note that the curvature Z̈(y) ∝ −h₁ cos ỹ − 4h₂ cos 2ỹ already shows a marked double well structure but the point representing the Li⁺ rows does not reach zero (h₂ < h₁/4) so that there is no inversion of the curvature. In other words, the observed central peak is a precursor of the genuine second classical rainbow peak.⁶⁰

For He where the same situation takes place in a semi-quantum regime, supernumerary rainbows are also visible in Fig. 16 of ref. 34, probably showing more details. With He along KCl[110], the ratio between ionic radii and lattice parameter is such that the same transition takes place close to the quantum regime, between 30 meV and 500 meV producing an apparent disorder in the diffraction charts reporting the evolution of the diffracted intensities^61–63 before a clear second “inner” rainbow structure becomes visible at E_⊥ = 0.9 eV.⁵³

It should be recalled that in the present paper, all decoherence and inelastic effects are embedded in the concept of effective line-profile measured here along the [Rnd] direction. In this description, the observed contrast weakening observed, for instance in Fig. 2 and 3 is reproduced by a fully coherent HCW model but with a line-profile becoming broader than the Bragg angle. We believe that this increasingly smooth profile should merge progressively into a fully classical profile in spite of a fully coherent model. As a result, the derivation of the corrugation amplitude should not be affected by the decoherence as it was shown when the analysis is based on quasi-elastic intensities.⁶⁴

7.2 Low energy part, attractive forces

For all projectiles and for both [100] and [110] directions, all data indicate a rapid increase of the rainbow angle Θ_rb when the effective energy E_⊥ approaches zero. The effect is well-known and, within the HCW, was modeled as the Beeby correction.⁶⁵ It assumes that the attractive forces are important at comparatively large distance from the surface. By replacing the initial wave vector k_⊥, associated with E_⊥, by , associated with , takes into account refraction of the rainbow angle following Snell's law. This also modifies the velocity and the phase shift k′z_c associated with a given path difference z_c (see e.g. the review⁶⁶).

Along the [110] direction we have focused on the rainbow angle position which is large enough to be pointed by hand and we have adjusted the value Θ_surf by hand so that the refracted curve passes through the data in Fig. 7 for Ar and Fig. 8 for Kr and Xe. Along the [100] the rainbow angle is comparatively small and we have used the very good agreement observed with the Bessel function I_m = J_m²(ζ) in Fig. 2 and 3, to derive precise values of the rainbow angle describing accurately the observed scattering distributions. For Ar, to recover the emission angle Θ_surf, we have applied the correction directly to the measured values in Fig. 7 and extracted the value Θ_surf. For Kr and Xe where the energy range is larger, we have adjusted the linear form by hand so that the refracted curve passes through the data in Fig. 11a and b respectively. All the values derived along both [100] and [110] directions are reported in Table 2 together with calculated values⁵⁵ of the well-depth D. This shows that the refraction effect using these values of the well-depth already provide a qualitative description of our data along both the [100] and [110] directions.

7.3 Model derivation of the PEL along [100]

The previous paragraphs have shown that the Beeby correction does significantly modify the apparent corrugation amplitude derived from the HCW. We have also shown that the same well-depth D of the mean planar potential can be applied along both directions⁶⁷ to produce a more realistic energy behavior. However, this latter was in part forced to follow a given arbitrary analytic form . This can also be related to an intrinsic weakness of the HCW, each measurement at an energy E_⊥ produces a value of ζ without any connection between them, i.e. ignoring the consistency between the iso-energy curves originating essentially from the quasi-exponential decay range of the PEL along the z direction.

A better consistency can be achieved by fitting all the data simultaneously with a model PEL so that the corrugation function z_E(y,z) are forced to be the iso-energy curves. Following ref. 24, we use a PEL based on binary potentials developed on a set of screened coulombic forms V(r) = Σ_i(A_i/r)e^−B_ir attached to atomic positions. Such forms are widely used in atomic collisions to describe the repulsive part.⁶⁸ For the attractive part, screened coulomb forms are also well-suited to describe the projectile polarization by the exponentially decaying Madelung field above the surface while van der Waals or Casimir-Polder terms are better modeled by Lennard-Jones term −B/r⁶ summing up to planar −B/4z⁴ dependence ending as ∝z⁻³ when retardation effects are taken into account.⁶⁹ For neon projectiles, the inclusion of the ∝z⁻³ parameter does affect the screened-coulombic attractive term but did not improve the quality of the fit²⁴ to the data which appeared to be sensitive mainly to the resulting well-depth D. We thus discard this planar term to keep the number of adjustable parameters as low as possible.

We will compare only the scattering data along the [100] direction where the data reduction to a single number ζ was established in Section 3 and illustrated in Fig. 2 and 3. The ζ(E_⊥) values derived from similar Bessel fits are displayed in Fig. 10 and 11 over a broad energy range (,). Note that the value ζ, Θ_rb or z_c, connected by eqn (3) represent the observed angular distribution so that, if present, refraction is embedded in these data.

Only four parameters are needed to describe the interaction potential along this direction, A₁,B₁ and A₂,B₂ with A₁ > 0 and A₂ < 0. The corrugation amplitude is calculated from the PEL averaged along [100] as z_c(E_⊥) = z_top − z_bot by solving numerically the equation PEL(y,z) = E_⊥ at the locations y = y_top and y = y_bot (see inset in Fig. 13). Using eqn (3) to connect z_c to Θ_surf and the value of D derived from the analytic mean planar form of the PEL to transform Θ_surf into the rainbow angle Θ_rb using Snell's law. This value is then compared with the measured ones reported in blue in Fig. 10 and 11a, b).

The corrugation amplitudes z_c (E_⊥) resulting from this least square fits are reported in Fig. 13 where published data for He³⁴ and Ne³⁶ have been fitted with the same procedure. No spurious increase is observed at low energy, instead a distinct leveling of the corrugation amplitude is present even when the contribution of the attractive forces on the projectile trajectories have been compensated.

This leveling is a direct consequence of the attractive forces, which have compressed all positive iso-energy lines below the one corresponding to E_⊥ = 0 rather than letting these lines expand smoothly to the vacuum. The popularity of the concept of corrugation amplitude is probably partly due to this flat behavior at low energy where most measurements have been performed. Fig. 13 also shows a significant decrease of the corrugation amplitude along the He–Xe sequence. In the present approach, the “size effect” is fully embedded in the repulsive binary potentials reported as A₁ and B₁ in Table 3. Neglecting attraction (i.e. forcing A₂ to zero) indeed pushes the turning point further away from the surface and geometrically reducing the corrugation amplitude as expected from a hard sphere model. This is not directly visible with our model PEL which indicates that the distance of the turning points does not change significantly. This is visible in Table 3 where z₀ represents the distance where the mean planar potential is zero. Since the attractive forces have a longer range, they translate the low energy iso-energy curves (E_⊥ ∼ D) closer to the surface preserving the reduced corrugation amplitude that would arise without attractive forces at E_⊥ ∼ D. The leveling region progressively merges a steady increase specific to the [100] direction where the region in between the atomic rows is empty. Note that the slope is sensitive to the attractive forces via the Rieder softness correction⁷⁰ rediscovered in GIFAD as a stiffness correction.⁶ Here again this is due to the compression of the equipotential lines by the attractive forces which increase the effective slope.

Table 3 Binary interaction parameters A_i, B_i adjusted to the data along the [100] direction (, in Fig. 10 and 11). The resulting mean planar potential, V(z) = Σ_iA_ie^−B_iz is close to a Morse potential with B₂ ≃ B₁/2. The values z₀ and z_e correspond to the V(z₀) = 0 and V(z_e) = D

	A 1	B 1	A 2	B 2	z 0	z e	D
He	127	1.71	−3.5	1.1	5.17	5.89	−10
Ne	537	1.96	−3.19	1.11	5.3	5.97	−10
Ar	340	1.64	−1.27	0.743	5.32	6.2	−55
Kr	1083	1.71	−2.26	0.765	5.62	6.52	−90
Xe	969	1.54	−3.0	0.679	5.78	6.74	−149

---

Returning to the analogy with the atomic force microscope, Fig. 13 allows a reformulation of the naive question “what is the size of the tip?” into better defined questions such as how many electrons are in the outer shell, what is the relative binding energy of the projectile relative to the surface conduction (or valence) band and what is the magnitude of the attractive forces.

It should be stressed that the fit should be comparatively robust, not only due to the redundancy provided by the reduced number of parameters but also because the fitted parameter ζ is a monotonous function where no oscillation is present avoiding the difficult situation where the fit is trapped in a local minima.

7.4 Perspectives

Extending our model to fit simultaneously the [110] data requires abandoning ζ which offers a reliable description of the data only along the [100] direction. The fitting parameter could be σ_m (or σ_ϕ) which is equivalent to ζ along [100] but is also perfectly defined along [110] and which can be extracted directly from the raw scattering profile. The data reduction procedure used above to extract the diffracted intensities I_m is a kind of de-convolution which becomes more and more difficult when the inelastic line width exceeds the Bragg angle. The reverse approach where the predicted intensities would be convoluted by the line profile recorded along the [Rnd] direction would be much more stable. Such a global fit of all directions would be even more robust and precise as it would be constrained by the supernumerary rainbow structure and could reveal the contribution of the Li⁺ ions. The short term perspective is to use the derived parameters A_i, B_i of the binary interaction potential as a seed for a fit by a quantum scattering code to produce values as exact as possible.

Over the longer term, fitting all noble gas data directly onto the surface electron density would be more instructive but requires an a priori knowledge of the van der Waals attraction and of the connexion between the surface electron density and the magnitude of the Pauli repulsion. The first point seems out of reach but the present study suggests that diffraction data are mainly sensitive to the well-depth D allowing a model description of attractive forces. Concerning repulsive contributions, a simple proportionality factor between the surface electron density and the magnitude of the Pauli repulsion was demonstrated for helium scattered from metals⁷¹ and is probably sound on LiF for the compact helium atom but becomes hardly justified for heavy atoms where the binding energies compare with the LiF surface work-function. Another interesting perspective would be to turn atomic diffraction sensitive to the chemical nature of the surface atoms. The semi-empirical Z. B. L.⁶⁸ or O. C. B.⁷² binary potential use similar screened coulomb forms with more parameters and are certainly a very good start to describe the interaction above a few eV but additional rules are needed to adapt the decay range of outer shell electrons to match, or not, that of the surface electron density.

8 Conclusion

In GIFAD the role of attractive forces was first overlooked,³⁰ in part because with He at low E_⊥, only a few diffraction orders are populated and the rainbow angle is not readily visible. The corresponding shift in the energy scale was first identified when comparing He diffraction on the β₂(2 × 4) reconstruction of the GaAs(001) surface with ab initio calculations.⁷³ The Beeby correction was used in GIFAD for He scattered from Graphene/SiC where the corrugation function was describing both the interatomic potential and the shape of the Moiré pattern.⁷⁴ Then the refraction effect was made model-free in Fig. 21 of ref. 34 by using the standard deviation σ_m or σ_Θ, bypassing the classical concept of rainbow angle irrelevant in the quantum regime where only one or two diffraction orders are populated. With neon on LiF[110] the reduced wavelength increases the number of diffraction orders making the refraction effect readily visible in the evolution of scattering profiles (Fig. 2 of ref. 24). By comparing with a quantum scattering code, the Snell's correction was shown to be a decent description of the refraction effect as already demonstrated for helium at low energy.⁶⁶ In the present paper, the refraction of Ar, Kr and Xe is shown to be quite large both along the [110] and [100] directions and is described consistently with a single value of the well-depth, in opposition with publications suggesting that van der Waals attraction is important only along the [110] direction.⁷⁵

The main interest of this semi-classical approach is that it is conceptually and numerically simple and the fit to all data takes place in a few seconds on an office PC. It could, for instance be used as a front-end to initiate the fit by quantum scattering codes offering both accuracy and exactitude.

In the specific case of a cosine corrugation function where the diffracted intensities are well-fitted by Bessel function I_m = J_m²(ζ), the parameter ζ is most often interpreted in terms of surface corrugation z_c = ζ/k_⊥ while, by definition it accurately describes the azimuthal profile P(m) = I_m and therefore its scattering width .⁴⁴ When trying to connect this value of ζ representing the observed scattering profile to the surface corrugation, the refraction by the attractive well has to be taken into account in the form of Snell's law in eqn (1) or with the Beeby correction: .

When refraction is taken into account, our results in Fig. 13 show that the corrugation amplitude becomes smaller when probed with heavier noble gas due to the increased magnitude of the Pauli repulsion, as if the turning point would be pushed away from the surface plane. This is not directly observed because the reduced “geometric” contribution is brought back towards the surface by the attractive forces in the form equivalent to a translation preserving the reduced corrugation and producing a distinct leveling at low energy E_⊥ ∼ D.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data will be made available upon any reasonable request.

Acknowledgements

We are grateful to Hynd Remita for the irradiation of the LiF samples by γ rays from the Cobalt source of the Institut de Chimie Physique at Orsay, favoring cleaving with large terraces. We are indebted to J. R. Manson for teaching us part of his vast background on atomic diffraction at thermal energies. We also acknowledge the support of Christian Brand for the stay of C. K. at Orsay and the enthusiastic participation of Alex Le Guen and Jaafar Najafi Rad during their internship in the group and the careful reading by A. J. Mayne. This work received support from LabEx PALM (ANR-10-LABX-0039-PALM) and Chinese Scholarship Council (CSC) Grant No. 201806180025.

Notes and references

1. A. Zugarramurdi and A. G. Borisov, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 86, 062903.
2. E. Pollak, P. Roncin, W. Allison and S. Miret-Artes, Phys. Chem. Chem. Phys., 2024, 26, 25501.
3. P. Roncin and M. Debiossac, Phys. Rev. B, 2017, 96, 035415.
4. P. Rousseau, H. Khemliche, N. Bundaleski, P. Soulisse, A. Momeni and P. Roncin, J. Phys.: Conf. Ser., 2008, 133, 012013.
5. J. R. Manson, H. Khemliche and P. Roncin, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 155408.
6. P. Pan, M. Debiossac and P. Roncin, Phys. Rev. B, 2021, 104, 165415.
7. P. Rousseau, H. Khemliche, A. G. Borisov and P. Roncin, Phys. Rev. Lett., 2007, 98, 016104.
8. A. Schüller, S. Wethekam, D. Blauth, H. Winter, F. Aigner, N. Simonovic, B. Solleder, J. Burgdörfer and L. Wirtz, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 82, 062902.
9. A. Muzas, F. Gatti, F. Martín and C. Díaz, Nucl. Instrum. Methods Phys. Res., Sect. B, 2016, 382, 49–53.
10. A. Zugarramurdi, M. Debiossac, P. Lunca-Popa, L. Alarcón, A. Momeni, H. Khemliche, P. Roncin and A. G. Borisov, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 012904.
11. A. P. Jardine, S. Dworski, P. Fouquet, G. Alexandrowicz, D. J. Riley, G. Y. Lee, J. Ellis and W. Allison, Science, 2004, 304, 1790–1793.
12. M. Debiossac, A. Zugarramurdi, P. Lunca-Popa, A. Momeni, H. Khemliche, A. G. Borisov and P. Roncin, Phys. Rev. Lett., 2014, 112, 023203.
13. B. S. Zhao, S. Schulz, S. Meek, G. Meijer and W. Schöllkopf, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 78, 010902.
14. R. I. Masel, R. P. Merrill and W. H. Miller, Phys. Rev. B, 1975, 12, 5545–5551.
15. B. Holst, G. Alexandrowicz and N. Avidor, et al. , Phys. Chem. Chem. Phys., 2021, 23, 7653–7672.
16. A. LeGrand and E. Greene, J. Chem. Phys., 1986, 84, 6483–6487.
17. J. M. Moix and E. Pollak, J. Chem. Phys., 2011, 134, 16.
18. M. S. Gravielle, J. E. Miraglia, A. Schüller and H. Winter, Nucl. Instrum. Methods Phys. Res., Sect. B, 2013, 317, 77–82.
19. J. E. Miraglia and M. S. Gravielle, Phys. Rev. A, 2017, 95, 022710.
20. L. Mattera, M. Rocca, C. Salvo, S. Terreni, F. Tommasini and U. Valbusa, Surf. Sci., 1983, 124, 571–582.
21. K. Rieder and W. Stocker, Phys. Rev. Lett., 1984, 52, 352.
22. E. Semerad, P. Sequard-Base and E. Hörl, Surf. Sci., 1987, 189, 975–983.
23. M. Minniti, C. Daz, J. F. Cuñado, A. Politano, D. Maccariello, F. Martn, D. Faras and R. Miranda, J. Phys.: Condens. Matter, 2012, 24, 354002.
24. M. Debiossac, P. Roncin and A. Borisov, J. Phys. Chem. Lett., 2020, 11, 4564–4569.
25. A. G. Stoll, D. L. Smith and R. P. Merrill, J. Chem. Phys., 1971, 54, 163–169.
26. R. Lahaye, A. Kleyn, S. Stolte and S. Holloway, Surf. Sci., 1995, 338, 169–182.
27. R. Lahaye, S. Stolte, S. Holloway and A. Kleyn, Surf. Sci., 1996, 363, 91–99.
28. B. Berenbak, S. Zboray, B. Riedmuller, D. Papageorgopoulos, S. Stolte and A. Kleyn, Phys. Chem. Chem. Phys., 2002, 4, 68–74.
29. D. Faras, C. Daz, P. Nieto, A. Salin and F. Martn, Chem. Phys. Lett., 2004, 390, 250–255.
30. H. Winter and A. Schüller, Prog. Surf. Sci., 2011, 86, 169–221.
31. C. Lim, J. C. Tully, A. Amirav, P. Trevor and M. J. Cardillo, J. Chem. Phys., 1987, 87, 1808–1816.
32. K. Gibson, S. Sibener, H. P. Upadhyaya, A. L. Brunsvold, J. Zhang, T. K. Minton and D. Troya, J. Chem. Phys., 2008, 128.
33. Y. Watanabe, H. Yamaguchi, M. Hashinokuchi, K. Sawabe, S. Maruyama, Y. Matsumoto and K. Shobatake, Eur. Phys. J. D, 2006, 38, 103–109.
34. M. Debiossac, P. Pan and P. Roncin, Phys. Chem. Chem. Phys., 2021, 23, 7615–7636.
35. M. Gravielle, A. Schüller, H. Winter and J. Miraglia, Nucl. Instrum. Methods Phys. Res., Sect. B, 2011, 269, 1208–1211.
36. M. Debiossac, P. Pan and P. Roncin, Phys. Chem. Chem. Phys., 2023, 25, 30966–30974.
37. P. Pan, J. N. Rad and P. Roncin, Rev. Sci. Instrum., 2022, 93, 1.
38. M. Debiossac and P. Roncin, Nucl. Instrum. Methods Phys. Res., Sect. B, 2016, 382, 36–41.
39. S. Lupone, P. Soulisse and P. Roncin, Nucl. Instrum. Methods Phys. Res., Sect. B, 2018, 427, 95.
40. P. Pan, M. Debiossac and P. Roncin, Phys. Chem. Chem. Phys., 2022, 24, 12319–12328.
41. M. Debiossac and P. Roncin, Nucl. Instrum. Methods Phys. Res., Sect. B, 2016, 382, 36–41.
42. J. Seifert, J. Lienemann, A. Schüller and H. Winter, Nucl. Instrum. Methods Phys. Res., Sect. B, 2015, 350, 99–105.
43. P. Pan, C. Kanitz, M. Debiossac, A. Le-Guen, J. N. Rad and P. Roncin, Phys. Rev. B, 2023, 108, 035413.
44. For instance, the width parameter for 4 keV Kr atoms was measured between 0.5° and 1.5° and modeled with a polynolial fit; σ_ϕ (mdeg) = 8 + 30 θ_i + 14 θ_i² (θ in deg). For 4 keV Xe, we convoluted the LG profile in eqn (2) by the beam profile (σ = 9 mdeg.) σ_ϕ = −1.5 + 46 θ_i + 21 θ_i².
45. Stricto sensu, there is only one specular spot corresponding to θ_f = −θ_i and ϕ_f = ϕ_i, we call here specular circle the circle corresponding to |θ_f| = |θ_i|.
46. T. Young, The Bakerian Lecture:: Experiments and Calculations Relative to Physical Optics, 2016.
47. P. Laven, J. Quant. Spectrosc. Radiat. Transfer, 2004, 89, 257–269.
48. H. Chow and E. Thompson, Surf. Sci., 1976, 59, 225.
49. S. Miret-Artes and E. Pollak, J. Phys. Chem. Lett., 2017, 8, 1009–1013.
50. U. Garibaldi, A. Levi, R. Spadacini and G. Tommei, Surf. Sci., 1975, 48, 649–675.
51. A. Schüller and H. Winter, Phys. Rev. Lett., 2008, 100, 097602.
52. A. Schüller, S. Wethekam, A. Mertens, K. Maass, H. Winter and K. Gärtner, Nucl. Instrum. Methods Phys. Res., Sect. B, 2005, 230, 172–177.
53. U. Specht, M. Busch, J. Seifert, A. Schüller, H. Winter, K. Gärtner, R. Wodarczyk, M. Sierka and J. Sauer, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 125440.
54. Note that for a cosine corrugation function where I_m = J²_m(k_⊥z_c), σ²_m = Σm²I_m is analytic: ².
55. H. Hoinkes, Rev. Mod. Phys., 1980, 52, 933.
56. E. Pollak and S. Miret-Artés, J. Phys. Chem. C, 2015, 119, 14532–14541.
57. W. Allison, S. Miret-Artés and E. Pollak, Phys. Chem. Chem. Phys., 2022, 24, 15851–15859.
58. As soon as the beam direction is not aligned with the crystal axis, perturbation theory56 clearly outlines the crucial role of the decay-range of the electron density also responsible for softwall corrections.
59. A. Momeni, P. Soulisse, P. Rousseau, H. Khemliche and P. Roncin, e-J. Surf. Sci. Nanotechnol., 2010, 8, 101–104.
60. S. Miret-Artés and E. Pollak, Surf. Sci. Rep., 2012, 67, 161–200.
61. E. Meyer, Habilitattionsschrift PhD dissertation, Humboldt Universität Berlin, Germany, 2015.
62. G. A. Bocan and M. S. Gravielle, Nucl. Instrum. Methods Phys. Res., Sect. B, 2018, 421, 1–6.
63. M. del Cueto, A. Muzas, F. Martn and C. Daz, Nucl. Instrum. Methods Phys. Res., Sect. B, 2020, 476, 1–9.
64. G. Anemone, A. Al Taleb, W. Hayes, J. Manson and D. Faras, J. Phys. Chem. C, 2017, 121, 22815.
65. J. L. Beeby, J. Phys. Chem. C, 1971, 4, L359.
66. D. Farías and K. Rieder, Rep. Prog. Phys., 1998, 61, 1575.
67. When comparing calculations or models to diffraction, the attractive part D is sometime used as an adjusting parameter but then, D should be the same along all directions. The conclusion that the VdW forces are important⁷⁶ along [110] but not along [110]⁷⁵ is surprising. These works contain no reference to refraction effect^24,66,77 or to the Beeby correction⁶⁵.
68. J. F. Ziegler and J. P. Biersack, Treatise on heavy-ion science: volume 6: astrophysics, chemistry, and condensed matter, Springer, 1985, pp. 93–129.
69. J. Lecoffre, A. Hadi, M. Bruneau, C. Garcion, N. Fabre, E. Charron, N. Gaaloul, G. Dutier and Q. Bouton, Phys. Rev. Res., 2025, 7, 013232.
70. K. H. Rieder and N. Garcia, Phys. Rev. Lett., 1982, 49, 43–46.
71. M. Manninen, J. K. Nørskov, M. J. Puska and C. Umrigar, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 2314–2316.
72. D. O’connor and J. Biersack, Nucl. Instrum. Methods Phys. Res., Sect. B, 1986, 15, 14–19.
73. M. Debiossac, A. Zugarramurdi, H. Khemliche, P. Roncin, A. G. Borisov, A. Momeni, P. Atkinson, M. Eddrief, F. Finocchi and V. H. Etgens, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 155308.
74. M. Debiossac, A. Zugarramurdi, Z. Mu, P. Lunca-Popa, A. J. Mayne and P. Roncin, Phys. Rev. B, 2016, 94, 205403.
75. G. A. Bocan, H. Breiss, S. Szilasi, A. Momeni, E. Staicu Casagrande, E. A. Sanchez, M. S. Gravielle and H. Khemliche, Phys. Rev. B, 2021, 104, 235401.
76. G. A. Bocan, H. Breiss, S. Szilasi, A. Momeni, E. S. Casagrande, M. S. Gravielle, E. A. Sanchez and H. Khemliche, Phys. Rev. Lett., 2020, 125, 096101.
77. J. G. Mantovani and J. R. Manson, Surf. Sci., 1982, 120, L487–L490.

---