Atom-Diffraction from Surfaces with Defects: A Fermatian, Newtonian and Bohmian Joint View

Bohmian mechanics, widely known within the field of the quantum foundations, has been a quite useful resource for computational and interpretive purposes in a wide variety of practical problems. Here, it is used to establish a comparative analysis at different levels of approximation in the problem of the diffraction of helium atoms from a substrate consisting of a defect with axial symmetry on top of a flat surface. The motivation behind this work is to determine which aspects of one level survive in the next level of refinement and, therefore, to get a better idea of what we usually denote as quantum-classical correspondence. To this end, first a quantum treatment of the problem is performed with both an approximated hard-wall model and then with a realistic interaction potential model. The interpretation and explanation of the features displayed by the corresponding diffraction intensity patterns is then revisited with a series of trajectory-based approaches: Fermatian trajectories (optical rays), Newtonian trajectories and Bohmian trajectories. As it is seen, while Fermatian and Newtonian trajectories show some similarities, Bohmian trajectories behave quite differently due to their implicit non-classicality.


Introduction
In the last several decades, there has been a fruitful and beneficial transfer of the ideas involved in David Bohm's formulation of quantum mechanics [1][2][3][4] from the domain of the quantum foundations to the arena of the applications [5][6][7][8][9]. The conceptual and mathematical background provided by Bohmian mechanics [3,4,[10][11][12] has become a resourceful tool to investigate quantum problems from an alternative viewpoint regardless of the always ongoing hidden-variable debate. This includes both fresh interpretations to (known and also new) quantum phenomena and novel implementations of alternative numerical algorithms to tackle them [13]. The essential ingredients of Bohmian mechanics have also inspired methodologies and descriptions aimed at providing effective trajectory-like explanations of wave phenomena beyond the quantum realm [14]. For instance, Bohmian-like trajectories have been synthesized from experimental data with light [15], while Bohmian-like behaviors have been recreated in classical fluid-dynamics experiments [16][17][18][19][20][21][22].
Getting back to quantum mechanics, one of the advantages of Bohmian mechanics is, perhaps, its capability to put on the same level quantum and classical analyses or descriptions of the same physical phenomenon by virtue of the concept of trajectory, well defined in both contexts. Now, because Bohmian trajectories are in compliance with quantum mechanics, they can be considered to be at a descriptive level above classical trajectories. Thus, an interesting question that naturally arises is how much or what kind of information is kept when passing from a descriptive level to the next one. This is precisely the question addressed here. To this end, a realistic working system is considered, although it is still simple enough to provide a clear answer to the question. Specifically, the phenomenon analyzed is the helium-atom diffraction from a carbon monoxide molecule (CO) adsorbed on a platinum (111) single-crystal surface. This is a system that has been extensively studied in the literature both experimentally and theoretically [23][24][25][26][27][28][29][30][31], even from a Bohmian viewpoint [32,33]. An appropriate description and knowledge of the CO-Pt(111) interaction is important to the understanding of the role of Pt as a catalyst of the electrochemical oxidation of the CO, with industrial and technological applications (of course, an extensive literature on other analogous systems is also available). In order to determine such an interaction, one of the experimental techniques employed is the He-atom diffraction (or scattering) at low energies (typically energies are between 10 meV and 200 meV, about three orders of magnitude below the range for low-energy electron diffraction). This is a rather convenient tool to investigate and characterize surfaces at relatively low energies with neutral probes, which provides valuable information about the surface electronic distribution without a damage of the crystal-the atoms remain a few Ångstroms above the surface, in contrast to low-energy electron diffraction, where electrons penetrate a few crystal layers, strongly interacting with the crystal atoms.
As is well-known, when a matter wave is diffracted by a crystal lattice, either by reflection (the case of He atoms or low-energy electrons) or by transmission (the case of neutron diffraction or high-energy electrons), the resulting spatial intensity distribution is characterized by a series of maxima along the so-called Bragg directions. The lattice structure can be determined from these characteristic patterns by means of a convenient modeling. Sometimes, however, these patterns have distortions due to a break of the translational symmetry typical of a perfect lattice. This symmetry-breaking can be induced by intrinsic thermal (lattice) atom vibrations (phonons) or by the presence of different types of defects randomly distributed across the lattice [34], for instance, adsorbed particles (adsorbates). In the case of a periodic surface, the presence of an adsorbate on top of it translates into a blurring of the well-defined Bragg features and the appearance of broad intensity features. This diffuse scattering effect [24] is analogous to the image distortion produced by a rough mirror-the larger the number of CO adsorbates on the clean Pt surface, the larger the distortion with respect to neat Bragg-like diffraction intensity peaks.
Instead of considering a large number of randomly distributed CO adsorbates, we are going to focus on the effects produced by a single isolated CO adsorbate. Moreover, we shall not focus on the quantitative description of the diffraction process itself, but on the analysis of how the features associated with a given descriptive level manifest on the next level, which is assumed to be more refined and, therefore, accurate. Accordingly, we will see that although such features are transferred from the model characterizing one level to the model corresponding to the next level, it not always easy to establish an unambiguous one-to-one correspondence. In simple terms, appealing to a biological metaphor, it is like considering a body and its skeleton. The skeleton is the structure upon which the body rests. However, although the body reveals some features of the skeleton (we can perceive the position of some bones under our flesh), it is a much more complex super-structure. In particular, here the problem considered is approached at three descriptive levels: • The Fermatian level, which refers to the analysis of the problem assuming a bare hard-wall-like (fully repulsive) model to describe the He-CO/Pt(111) interaction. Because the trajectories here are of the type of sudden impact (free propagation except at the impact point on the substrate wall, where the trajectory is deflected according to the usual law of reflection), they are going to be straight-like rays, as in optics (this is why it is referred to as Fermatian). • The Newtonian level, where the He-CO/Pt(111) interaction is modeled in terms of a potential energy surface that smoothly changes from point to point. This model has a repulsive wall that avoids He atoms to approach the substrate beyond a certain distance (for a given incidence energy), and an attractive tail that accounts for van der Waals long-range attraction. The existence of these two regions, repulsive and attractive, gives rise to an attractive channel around the CO adsorbate and that continuous along the flat Pt surface, inducing the possibility of temporary trapping for the He atoms.
• The Bohmian level is the upper one and, to some extent, makes an important difference with the previous models because here the trajectories are not only dependent on the interaction potential model, but also on the particular shape displayed by the wave function at each point of the configuration space at a given time (the "guiding" or "pilot" wave).
From a physical view point, we are going to focus on three different aspects or phenomena, namely reflection symmetry interference, surface rainbows and ısurface trapping. As it will be seen, all these aspects have a manifestation in the corresponding diffraction intensity patterns. Reflection symmetry interference is the mechanism considered to explain the oscillations displayed by the intensity pattern at large diffraction angles [24,26], and is based on the hypothesis that the diffracted wave can be assumed to be constituted by two interfering waves. In terms of trajectories, as we shall see, this means that there are pairs of homologous paths (either Fermatian or Newtonian), which, when using semiclassical arguments in terms of the optical concept of paths difference, explain the appearance of such type of interference [24,26]. Rainbow features arise as a consequence of the local changes in the curvature of the interaction potential model [27,35], which give rise to accumulations of classical trajectories along some privileged directions (rainbow deflection angles), but that quantum-mechanically leave some uncertainties when we look at the corresponding diffraction patterns [29,32]. Finally, surface trapping along the clean Pt surface arises for some deflections from the adsorbate at grazing angles. This phenomenon takes place when, by virtue of the interaction with the adsorbate, the He atoms lose too much energy along their perpendicular direction, transferring it to their parallel degree of freedom (perpendicular and parallel are defined with respect to the clean Pt surface). The energy associated with their perpendicular degree of freedom becomes negative, while the parallel energy gets larger than the incident one (by conservation), which ends up with the atom moving in the form of jumps along the surface until it encounters another adsorbate that, by means of the reverse process, can be used to release the atom from the surface [36].
According to the above discussion, this work has been organized as follows. Details about the interaction potential models considered to determine the different trajectory dynamics are provided and discussed in Section 2. In this section, a brief discussion is also given concerning the numerical approaches used to determine the calculations shown here. The diffraction intensity patterns for both the hard-wall model and the realistic interaction potential are shown and discussed in Section 3. Section 4 is devoted to the description and analysis of the different types of trajectories (Fermatian, Newtonian and Bohmian), emphasizing the features that are both common to all approaches and also their main differences. To conclude, a series of final remarks are summarized in Section 5.

Potential Model and Computational Details
Interaction potential models for systems like the one we are dealing with here are determined from information extracted from experimental diffraction patterns. Thus, let us consider that, as it is typically done, the intensity distribution in these patterns is specified in terms of the transfer of He-atom momentum from the incidence direction to the direction parallel to the surface or, in brief, parallel momentum transfer [24], i.e., where θ d and θ i are the diffraction and the incident angles, respectively, and k i = √ 2mE i /h is the incident wavenumber. Then, diffraction features in the large-angle region of the pattern (large parallel momentum transfers) typically convey information about the repulsive part of the interaction, while the attractive part has a more prominent influence on the small-angle region (low parallel momentum transfers). From these potential models, it is possible to estimate the effective size of the adsorbed particles [23,24] as well as the local curvature of the surface electron density, which additionally may induce the presence and contribution of rainbow-like features [35,37] whenever the interaction potential consists of a short range repulsive region followed by a longer range attractive one (accounting for the van der Waals interaction between the adsorbate-surface system and the impinging neutral atom).
The model considered here, proposed by Yinnon et al. [27], gathers the above mentioned features and nicely fits the experimental data [27][28][29]. This potential model consists of two terms: where V Pt(111) (z) and V CO (r) describe, respectively, the separate interaction of the He atom with the clean Pt surface and the adsorbate. The functional form of these two potential functions are a Morse potential for the He-Pt(111), and a Lennard-Jones potential for the He-CO, with r = x 2 + y 2 + z 2 . The reference system for the potential (2) is centered at the CO center of mass, with D = 4.0 meV, α = 1.13 Å −1 , z m = 1.22 Å and = 2.37 meV [28,29]. With this, z accounts for the position of the He atom along the perpendicular direction with respect to the clean Pt surface, while x and y describe its parallel motion. Because of the rotation symmetry around the axis x = y = 0 exhibited by the interaction potential (2), for an illustration, in the contour plot displayed in Figure 1a only one half of the transverse section along the plane y = 0 is shown. In it, negative and positive energy contours are denoted, respectively, with blue solid line and red dashed line (due to symmetry, only a half of the potential is shown). On the right-hand side, panels (b)-(d) show different transverse sections of the potential to better appreciate the effect of the local curvature along the z direction at three different distances from x = 0 (for x = 0 Å, x = 3.31 Å and x = 6.35 Å, respectively), which give an idea of the well-depth on top of the adsorbate (about 2.96 meV), at the intersection of the adsorbate with the flat Pt surface (6.37 meV), and on top of the flat surface far from adsorbate (4 meV), respectively. The existence of an attractive well around the adsorbate and also along the surface is going to induce temporary trapping both classically and quantum-mechanically-only the presence of another neighboring adsorbate may help to remove such trapping. Since far from the influence region of the adsorbate the He-Pt(111) interaction only depends on the z coordinate, the trapped motion will be ruled by the well of the Morse function, being free along the x direction. The resulting motion is thus a combination of jumps along the z direction with a uniform motion parallel to the Pt surface, with an average speed larger than √ 2E i /m, with E i being the incident energy (notice that the energy along the x direction has to be larger than along the z direction in order to achieve negative values along the latter and, hence, trapping). Classically, if the energy along the z direction is given by where E z = E i − E x < 0 and m is the He-atom mass, it is easy to show that the turning points will be located at This motion is anharmonic, with frequency The length of each jump along the x direction can be easily estimated from Equation (7) according to the relation where p x = mv x and τ = 2π/ω. As can be noticed, at the threshold E x = E i , the jump length becomes infinity, i.e., the trajectories leave the adsorbate being and remaining parallel to the clean Pt surface. The energy difference between consecutive repulsive/attractive contour levels (red dashed lines/blue solid lines) is 10 meV/1 meV. The thick black solid line denotes the repulsive boundary for an approximate hard-wall model set for an incidence energy of 10 meV (see text for details). In the right-hand side panels, energy profiles along the z direction for: (b) x = 0 Å; (c) x = 3.31 Å and (d) x = 6.35 Å. In these panels, the total interaction potential is denoted with black solid line, while red dashed and blue dash-dotted lines refer to the Morse and Lennard-Jones contributions, respectively.
Quantum-mechanically, the trapped portions of the wave packet will somehow contain information about the bound states of the Morse function, with eigenenergies given by [38] with being the harmonic frequency resulting from approximating the Morse potential to a harmonic oscillator. From the condition to determine the total number of bound states, i.e., ∆E n =n+1,n = E n − E n ≥ 0, it is found that, for the parameters here considered, there are only three bound states, namely the ground state plus two excited ones: E 0 = −2.53 meV, E 1 = −0.60 meV and E 2 = −3 × 10 −3 meV. If we assume E z = E n , then we obtain ∆ (0) x ≈ 23 Å and ∆ (2) x ≈ 320 Å, respectively, according to Equation (8).
The potential model (2) will be used to investigate the behavior of Newtonian and Bohmian trajectories when the He atoms are influenced by it. For the Fermatian approach, we shall consider a rather crude approximation to this potential, which consists of substituting it by a purely repulsive infinite barrier, which will be referred to as the (repulsive) hard-wall model. This barrier is an approximation to the equipotential V(r) = 10 meV (black thick solid line in Figure 1a), so that V HW (r) = ∞ for any position r such that V(r) ≥ 10 meV and zero everywhere else (i.e, for all r such that V(r) < 10 meV). An optimal fit to this equipotential renders a radius a = 2.86 Å for the adsorbate and a position z r = 0.28 Å of the clean Pt surface above the CO center of mass. Since the Bohmian treatment of this model is not included here, it is worth mentioning a recent analysis of an analogous system by Dubertrand et al. [39], where the authors consider the Bohmian description of quantum diffraction by a half-flat surface (simplified by a half-line barrier), earlier considered by Prosser [40] in the context of the solution of Maxwell's equations to the problem of diffraction by a semi-infinite conducting sheet [41].
One of the advantages of the model introduced in this section is that it can describe the presence of a localized adsorbed particle on a surface (which can be regarded as a point-like defect [24]), but also a row of adsorbates aligned, for instance, along the y-axis (a linear-like defect [23]). In the first case, we have radial symmetry with respect to the z-axis, while the latter is characterized by axial symmetry, along the y-axis. From the viewpoint of a trajectory-based description, there is no difference between one case and the other, since what happens in one half of a transverse section also happens in any other section (which can be reconstructed either by rotation symmetry around the z-axis or by translational symmetry along the y-axis and/or mirror symmetry with respect to the yz plane). The difference between both models relies on the way how the trajectories distribute spatially and therefore how many of them lay within a certain solid angle, independently of whether the trajectories are classical (Fermatian or Newtonian) or quantum-mechanical (Bohmian). Having said this, since we are interested in comparing the behavior of different types of trajectories, from now on, the discussion will turn around the two-dimensional description, which is analogous to consider an axial-symmetric system.
Regarding the conventions used here [42], the incidence angle, θ i , for Fermatian and Newtonian trajectories is defined as the angle subtended between the incident direction of a given trajectory and the normal to the clean Pt surface. Dynamically speaking, this translates into an effective way of how much momentum is provided initially to each direction, i.e., with p i =hk i = √ 2mE i . In particular, in the calculations presented and discussed below, we have considered two values of the incident energy, namely E i = 10 meV and 40 meV. In terms of the de Broglie wavelength, λ dB = h/ √ 2mE i , with h being Planck's constant, these energies correspond to λ dB = 1.43 Å and 0.72 Å, respectively. The relations (11) are also used for the quantum analysis, where the incident wave function is launched with a momentum in compliance with these expressions. The deflection (or outgoing) angle, θ d , on the other hand, is defined as the angle subtended by the normal and the deflection direction for the corresponding trajectory-in the case of wave-function descriptions, this angle is going to be denoted as the diffraction angle -, although an analogous definition in terms of the momentum with which the particle is deflected can also be used. Once the incidence angle is established, depending on the initial position assigned to the trajectories, they will behave in a way or another. To characterize the trajectories according to their initial positions, it is common to refer them to the impact parameter, which in the present context is defined as the impact position on the clean surface in absence of interaction. In periodic surfaces, the range for impact parameters (b) is typically established in terms of the lattice parameter (the unit cell length) [42,43].
Here, because the presence of the adsorbate breaks the periodicity of the clean surface, we need to redefine this range in an alternative and slightly different way. Specifically, this range is taken as a portion of surface that covers the extension of the adsorbate and goes well inside the region where the flat surface potential is already not influenced by the adsorbate attractive tail. With the potential function used here, this means is satisfied by impact parameters taken within the range [−10.6, 10.6] Å. Accordingly, the initial positions are specified as with b ∈ [−10.6, 10.6] Å, and where z 0 = 10.27 Å has been chosen far from the surface, such that V(r) ≈ 0 (the same holds for an incident wave function, whose probability density has to be far from the influence region of the adsorbate).
With the above definitions, the computation of Fermatian trajectories is trivial, since we only need to know their incidence direction and, from it, the point on the substrate (adsorbate or flat surface) where they will feel the impact. At such a point, applying the law of reflection, we readily obtain the deflection angle and therefore the deflected part of the trajectory. Fermatian trajectories are just a pure geometric issue and, as it will be seen in next section, the corresponding quantum calculations are just analytical, so they do not imply high computational demands. For Newtonian trajectories, however, the computational task is more refined, since the action of the interaction potential introduces important changes in the curvature of the trajectories when they are approaching the substrate. Nevertheless, the computational demand is still low, since such trajectories can be readily obtained by integrating Newton's equations (actually, Hamilton's equations) with a simple fourth-order Runge-Kutta algorithm using the above momentum and position values, Equations (11) and (12), respectively, as initial conditions.
The numerical computation of the wave-function evolution and the Bohmian trajectories is, however, more subtle, since it implies the solution of a partial differential equation. In this case, integration has been carried out by means of the second-order finite-difference algorithm [44], making use of the fast Fourier transform to compute the kinetic part of the operator [45]. For the initial wave function, we would like to simulate an incident nearly plane wave, which mimics a highly collimated He-atom beam. Numerically, we can recreate this situation by considering a quasi-monochromatic wave function or wave packet that covers the substrate well beyond the effective size of the adsorbate. This can be done by linearly superimposing a large number of Gaussian wave packets, which in our cases amounts to considering 250 Gaussian wave packets [46], where the spreading of each wave packet is 0.84 Å along the x direction and 2.65 Å along the z direction. With these conditions, the wave function reaches the surface with almost no increase of its size, which is launched from an average position along the z direction z 0 = z 0 = 10.27 Å and normal incidence conditions (again, for visual clarity). In order to ensure an optimal overlapping along the x direction, the centers of the wave packets are separated a distance of 0.21 Å. For simplicity both here and also with the hard-wall model, the quantum calculations have been performed at normal incidence conditions (θ i = 0 • ), although this does not diminish the generality of the results presented.
Regarding the computation of the Bohmian trajectories, they are synthesized on the fly from the wave function. That is, the wave function is made to evolve for a small time interval dt, and then the trajectories are propagated from their actual position to the new one with the phase information provided by the updated wave function. The equation of motion that rules this behavior is the guidance Equation (24) (see Section 4.3 for further details on this equation of motion). Since the value of the wave function and its derivatives is known only on the knots of the numerical grid, the guidance equation has to be solved with the aid of numerical interpolators, which render the values required at any other point other than a grid knot with a reliable accuracy. With these values, the equation of motion is solved by means of a Runge-Kutta algorithm, as in the case of Newtonian trajectories, although the degree of accuracy required is higher, particularly due to the appearance of nodal structures. As for the initial conditions, they have been chosen along lines parallel to the flat surface, at different constant distances from the latter and taking the value z 0 = z 0 as a reference. Specifically, three of these lines have been taken above this value [z(0) > z 0 ] and another three have been taken below [z(0) < z 0 ], i.e., closer to the substrate (see Section 4.3 for further details).

Diffraction from a Repulsive Hard-Wall Potential
In the case of the two-dimensional (axial-symmetric) version of the hard-wall model described in Section 2, the diffracted wave far from the adsorbate can be obtained from the exact (analytical) asymptotic solution to the problem of the diffraction from a cylinder [47], with r → ∞, and where is an outgoing wave vector pointing along an arbitrary diffraction direction-the two wave vectors are expressed in terms of their parallel and perpendicular components (k i,x and k d,x , and k i,z and k d,z , respectively) on purpose. In this expression, the first term is the contribution from the direct wave and the second term accounts for the diffraction caused by the defect itself (the minus sign in the first contribution arises from the fact that the incidence direction is considered to be negative).
If instead of a cylinder we have a half of it on top of a flat surface, the solution (13) has to include the effect of the reflection from the flat surface, i.e., which has to satisfy the hard-wall boundary conditions The first two terms in (14) satisfy this condition when z = 0, as expected on the flat surface. On the other hand, for the third term to satisfy these boundary conditions it is necessary that the diffraction (or scattering) amplitude, f , is given by two contributions, where the first term describes direct reflection from the adsorbate and the second, a double reflection from the adsorbate and then the flat surface. Both amplitudes can be recast in terms of the difference between the diffraction and incidence angles. In the first case, this can be readily seen; in the latter, a similar result is obtained after assuming collisions with a full cylinder, because then the diffraction angle is θ d = π − θ d . In addition, note that the symmetry displayed by Equation (15) for the specific case θ i = 0 is analogous to the antisymmetry condition arising in fermion-fermion collisions [48], where the symmetrized amplitude for two fermions with 1/2-spin in a triplet state has the functional form Analytically, in the short-wavelength limit (ka → ∞) and for a cylindrical defect, the diffraction amplitudes in Equation (15) are of the form [47], where the first term, known as the illuminated face term, accounts for the backward scattering of the wave and the second one describes the Fraunhofer diffraction. The final analytical expression for the diffraction amplitude f is obtained by considering the symmetrization condition (15) in these results.
The diffraction intensity for the axial-symmetric, fully repulsive hard-wall model is shown in Figure 2a as a function of the parallel momentum transfer (1), for E i = 10 meV and normal incidence (θ i = 0 • ). In the same figure, the intensities associated with the illuminated face term and the Fraunhofer diffraction are also displayed separately (red dotted line and blue dashed line, respectively) in order to get a better idea at a quantitative level of their respective contributions to the total pattern. As it can be seen, for small momentum transfers (small diffraction angles), the leading term is the Fraunhofer one, which decreases fast as the momentum transfer increases (as (∆K) −2 ). On the contrary, the illuminated-face term, together with its mirror image, becomes the leading contribution for large momentum transfers (large diffraction angles). The type of oscillations generated by these two terms, the illuminated-face one and its mirror image, give rise to the reflection symmetry phenomenon [24], which explains why the diffraction pattern does not decay for large momentum transfers, as happens, for instance, in simpler cases of diffraction by a wire or a slit. This behavior is observed regardless of the incident energy, as can be noticed from the intensities displayed in Figure 2b for E i = 40 meV (and also normal incidence).

Diffraction from the Potential Model (2)
The quantum treatment for the potential model (2) does not admit analytical solutions, as it is the case of the hard-wall model of Section 3.1, which provides us with an asymptotic analytical solution of the diffraction far from the interaction region. A way to tackle the problem is by using a numerical wave-packet propagation method, as described in Section 2, which renders a description of the diffraction phenomenon in real time, i.e., providing us with direct information on the time-evolution of the He-atom wave function, of particular interest in the region where the interaction between the He atom and the substrate is stronger. Hence, although we lose analyticity, we gain insight on the dynamical process in the interaction region.
Accordingly, the evolution of the probability density as it approaches the adsorbate and then gets diffracted is shown in Figure 3 at three different instants of its evolution for E i = 10 meV and normal incidence (taking advantage of the mirror symmetry with respect to x = 0 due to the normal incidence, only a half is plotted for simplicity). In panel (a), we observe the appearance of circular wavefronts (ripples) around the adsorbate due to the interference produced by the overlapping of the part of the wave function that is still approaching the adsorbate with the part that is already being diffracted. In panel (b), the whole of the wave function is interacting with the substrate (at about 1.5 ps). In this case, there are circular wavefronts around the adsorbate, as before, but also additional plane wavefronts produced by an analogous interference process associated with the portion of the flat Pt surface reached by the incident wave function. The superposition of the circular wavefronts with the planar ones generates around the adsorbate a web of maxima (see inset for a more detailed picture), which gets weaker and blurred as we move far from the adsorbate. The periodicity of the web of maxima is related to the incoming He-atom de Broglie wavelength, λ dB = 1.43 Å, although with some distortions due to the presence of the attractive region of the interaction. As the wave function further evolves, we can observe a superposition of two contributions, as seen in panel (c). One of them is the reflection of a nearly square function, which starts displaying the typical Fresnel or near-field features of such functions when they arise from single slit diffraction. This contribution is precisely the illuminated face term that we saw in Section 3.1. The other contribution, which distributes around the adsorbate, is going to be related to the Fraunhofer diffraction and also to other features, such as trapping of the wave along the attractive well near the Pt surface. Nonetheless, notice that the final diffraction pattern, as in the case of the hard-wall model (see Equation (17)), is a superposition of the two contributions (diffraction amplitudes), and not only the direct sum of their probabilities. On the left-hand side, contour plots illustrating three different instants of the evolution of the probability density near the surface for incidence conditions θ i = 0 • and E i = 10 meV: (a) when the density starts being influenced by the adsorbate; (b) when the density is totally interacting with the substrate (i.e., with both the adsorbate and the flat Pt surface) and (c) when the density starts leaving the substrate. In panel (d), on the right-hand side, plot of the probability density far from the influence of the adsorbate (t = 11 ps). Arrows and capital letters denote different diffraction directions to be identified in the intensity plot displayed in Figure 4a: A i : directions identifying interference features associated with the superposition of the circular and planar wavefronts, contributing to the central maxima of the intensity pattern; B i : associated with features arising from the reflection symmetry interference phenomenon; C: surface trapping.
Asymptotically, far from the influence of the classical interaction, as shown in Figure 3d, it is possible to observe traits related to the evolution of the three parts of the scattered wave, which can be eventually associated with the peaks characterizing the corresponding intensity pattern (see Figure 4a). Thus, interferences arising from the superposition of the circular and planar wavefronts give rise to the central features of the intensity pattern (denoted with directions labeled with A i in the plot). However, there are also other types of peaks, which propagate along the directions denoted as B i and, as will be seen below, can be associated with the reflection symmetry interference phenomenon. Both types of peaks, A i and B i , implicitly carry information about the classical rainbow (see Section 4.2) in a sort of global fashion, since this phenomenon does not manifest with a particular well-defined peak, in general. On the other hand, it is also possible to observe the presence of trapping (C), which is related to the lower part of the wave function that keeps moving close and parallel to the clean Pt surface. In order to quantify the effects produced by the diffraction process in the incoming wave, we represent the diffracted wave function as a linear combination of plane waves, where the elements S(k i,x , k d,x ) provide the probability amplitudes associated with the change of parallel momentum ∆K = k d,x − k i,x . These elements are determined by projecting the numerically computed diffracted wave onto this expression, from which the probability |S(k i,x , k d,x )| 2 is obtained to detect atoms that have exchanged a given amount ∆K of parallel momentum is obtained (i.e., the reflection coefficient). With periodic lattices, this calculation is typically performed within a single unit cell; here, because of the lack of periodicity, the calculation involves an artificial cell of about 53 Å, which covers a region large enough as to include both the adsorbate and a good portion of the flat surface that is not influenced by artifacts related to the adsorbate curvature-this is appropriate to capture isolated signatures of trapping (otherwise, there could be some contamination from the wave scattered in other directions). This procedure, however, has an inconvenience: the intensity pattern includes a rather high contribution from mirror reflection from the flat Pt surface. In order to remove it, the incident plane wave, ranging from −∞ to ∞, is decomposed as a linear superposition of two contributions, one contained within the integration range R (φ) and another outside it (χ), i.e., where the subscript i labels functions from a given basis set established on purpose to construct the wave function. Because the χ waves essentially describe diffraction from the clean Pt surface, the last term in this expression is going to be a contribution with the form of a δ-function along the incident direction. Thus neglecting this contribution, if the diffraction process from the surface is altered, we make a projection of the wave on the plane-wave basis set, and the S-matrix elements can be recast as where S R 0 (k i , k d ) is the matrix element corresponding to the scattering of the wave by a flat unit cell. The intensity to be compared with the experiment is thus obtained from the differential cross section or differential reflection coefficient, This intensity is displayed in Figure 4a as a function of the parallel momentum transferred (instead of the deflection angle, θ d ) after the wave function has evolved for 5.5 ps after the maximum approach to the surface. By this time, the action of the adsorbate interaction potential on the wave function is already negligible. Comparing the solid line with the dashed line, we notice the effect of removing the contribution from the plane wave. Accordingly, using (21) is analogous to "smoothing" the diffraction pattern, where the probability distribution does not show well defined maxima because of the presence of the plane-wave contribution. The clean oscillations that we observe (solid line) arise from the interference between the circular wave fronts coming from the adsorbate and the plane wavefronts coming from reflection from the flat Pt surface. This "interaction" is more prominent for small values of ∆K, this being the reason why the pattern, before removing the plane-wave contribution, displays fast oscillations. Such oscillations, however, get weaker and even meaningless as ∆K increases, because the circular wavefronts become less affected by their overlapping with the flat outgoing wavefronts. This is analogous to the behavior already observed with the hard-wall model, in Figure 2a: for small values of ∆K, the dominant contribution was the Fraunhofer one, while, for larger values of ∆K, the leading one was the illuminated-face contribution.
The intensity maxima in Figure 4a, though, do not totally correspond with those in Figure 2a, even if their number is the same. In particular, notice that some of these maxima (A 1 and A 3 ) display a kind of "wings". In order to elucidate their origin, the same calculation has been repeated using a repulsive model for the adsorbate, which consists of removing the attractive part of the Lennard-Jones function (4). The results for this model are displayed in panel (b). Comparing both models, we find that everything is essentially the same, except precisely for the presence of such "wings". This result is indeed close to the one displayed in Figure 2a for the hard-wall model, although the maxima are wider, which can be associated with the presence of an attractive well around the substrate. In order to determine whether it is an effect or not linked to the incidence energy, the same analysis was repeated for E i = 40 meV and normal incidence. From the classical calculations, we conclude that such "wings" have to be associated with the presence of rainbows [35] (see Section 4.2), since this phenomenon is linked to the local curvature of the interaction potential around the adsorbate. In this regard, the attractive well around the adsorbate plays a key role, since its removal makes rainbow features to disappear even though there is still an attractive region along the clean Pt surface. Actually, on a more quantitative level, notice that, for instance, for E i = 10 meV, the rainbow appears for ∆K ≈ 1.89 Å −1 . In Figure 4b, this value of ∆K is close to the maximum of the second lobe, which would lead to a distortion of this maximum and the adjacent ones (where the traits A 1 and A 3 appear). In the case of E i = 40 meV, there is a classical rainbow at ∆K ≈ 1.13 Å −1 , which corresponds to the maximum of the second lobe and, in consequence, this lobe essentially disappears in the attractive model. Thus, the quantum manifestation of classical rainbow features does not necessarily mean a contribution in terms of given maximum in the intensity pattern [28,29], but it can also be in terms of "global" phenomenon that affects the whole pattern [35,49]. This has actually been a controversial point in the literature when assigning the origin of the different diffraction maxima [27][28][29].
Finally, regardless of whether the interaction model considers an attractive or repulsive adsorbate, and also independently of the value of the incident energy, we observe that the last lobe of the intensity plots in Figure 4 gathers information of both grazing deflection and trapping. In panel (a), for instance, this maximum has been label as B 2 + C. According to Figure 3d, the maximum for B 2 should appear at ∆K ≈ 4.16 Å −1 . On the other hand, in the same figure, for C, we should have a maximum parallel transfer. However, the trapped probability is indeed oscillating inside the well while it moves along the x direction, which makes the corresponding ∆K value to fluctuate. Thus, the probability amplitudes related to B 2 and C will display some interference, which is precisely what we observe in the last lobe of all calculations presented in Figure 4.

Fermatian Level
Although it is a rather crude approximation, the hard-wall model is quite insightful because it allows for explaining and understanding on simple terms the reflection symmetry interference phenomenon [23,24,27,47,50,51] as well as the conditions leading to trapping. As in geometric optics, the key element is the interpretation of wave phenomena in terms of the phase difference arising from two different but equivalent paths, where by "equivalent" we mean that both leave the surface with the same deflection (outgoing) angle, even if their journeys close to the surface are quite different (actually, it is this difference that generates the phase difference). These geometric rays are what we call here Fermatian trajectories with the purpose to highlight such optical connotation and, as seen below, in the particular problem we are dealing with here, there are always homologous pairs of such trajectories. These pairs are formed by one trajectory that undergoes a single bounce from the substrate before getting deflected and another trajectory that undergoes two bounces (one with the adsorbate and another with the flat surface).
In Figure 5, there is a set of Fermatian trajectories, F , of particular interest: they are separatrices that determine the boundaries for ensembles of Fermatian trajectories that display a particular behavior, that is, all the Fermatian trajectories confined within two adjacent separatrices are going to exhibit an analogous behavior. In general, it can be noticed that, for a given incidence angle (here, θ i = 20 • ), trajectories may display either a single collision (regions denoted with light gray) or double collisions (blue and purple regions, denoted with A, B and C) depending on their impact parameter. The deflection can then be forward (trajectories represented with solid line) or backward (dashed line). Moreover, there are also regions of geometric shadow (red region, denoted with S) that cannot be reached by any trajectory, except for normal incidence (θ i = 0 • ). This region covers an area of length which depends on θ i , vanishing ( = 0) for θ i = 0 • and reaching its maximum extension ( = ∞) for parallel incidence (θ i = π/2). Each trajectory in Figure 5 carries a label, which helps to identify behavioral domains. We have that any trajectory impinging on the substrate either to the left of F α or to the right of F β will undergo forward deflection; any other trajectory confined in between will be deflected backwards. Notice that F α and F β are the only two trajectories that are deflected along the normal, constituting themselves a pair of homologous trajectories, with F α displaying double collision (first with the flat surface and then with the adsorbate) and F β displaying a single collision (with the adsorbate). As for the other separatrices:

•
Trajectories to the left of F 1 or to the right of F 7 only interact with the clean Pt surface and hence their deflection and incidence angles are equal. These trajectories, plus F 5 only contribute to mirror reflection from the flat surface, only contributing the intensity for ∆K = 0, since θ d = θ i -hence, this contribution will be more prominent as the range of impact parameters increases. • Any trajectory between F 1 and F α is deflected in an angle that goes from θ i to 0 • as the impact parameter increases. The same deflection angles are found for trajectories between F β and F 5 , although here the trend is that the angle increases from 0 • to θ i as b increases. Here, we have two sets of pairs of homologous trajectories: trajectories from the first set undergo double collisions (first with the flat surface and then with the adsorbate) and trajectories from the latter only have a single collision (with the flat surface). For any of these pairs, the angular distance between their impact points on the adsorbate surface is π/2 − θ i , as can be seen in Figure 6a.

•
There are also pairs of homologous trajectories with deflection angles between θ i and π/2. These are the trajectories confined between F 5 and F 6 , with single collisions (with the adsorbate), and between F 6 and F 7 , with double collisions (first with the adsorbate and then with the flat surface). This second set corresponds to trajectories impinging on the adsorbate within the sector B. In this case, the angular distance between impact points is not a constant, but depends on the deflection angle as π/2 − θ d . This distance gradually vanishes as both trajectories approach R 6 and is maximum when the trajectories coincide with the separatrices R 5 and R 7 . A representative set is depicted in Figure 6b. • Trajectories F 2 and F 4 are both deflected backwards along the incidence direction, i.e., θ d = −θ i . Accordingly, trajectories between F α and F 2 are deflected between 0 • and −θ i after undergoing double collisions (first with the flat surface and then with the adsorbate), while trajectories between F 4 and F β (to the right of F 4 ) undergo single collisions. The angular distance between impact points of homologous pairs of trajectories is now π/2 − θ d , although not all trajectories between F 4 and F β have a correspondent between F α and F 2 . This is because the flat surface intersects the adsorbate surface at a distance z = z r above its center of mass instead of at z = 0. Thus, instead of reaching a maximum deflection of −θ i , we have θ max d = −θ i + (sin) −1 (z r /a), which is the deflection for the trajectory F 2 . An illustrative pair of homologous trajectories of this kind is displayed in Figure 6c.

•
The trajectory F 3 separates the sets of homologous trajectories that are backward deflected, with the second collision taking place from the flat surface. One set is confined within trajectories F 2 and F 3 , with double collisions (first with the adsorbate and then with the flat surface), and the other set, with single collisions, is delimited by F 3 and F 4 (trajectories to the left of F 4 ). Unlike the previous set of backward-scattered homologous pairs, here all trajectories are paired, with the angular distance between their impact points being π/2 − θ d , as before. A representative pair is displayed in Figure 6d. In general, independently of whether the above pairs of homologous Fermatian trajectories describe situations of forward or backward scattering, and also regardless of the incidence angle, double collisions arise whenever the impact on the adsorbate surface takes place between the point where this surface intersects the flat Pt surface and the point at which any impinging trajectory is deflected parallel to such a flat surface (at a given incidence). These conditions determine the regions labeled with A, B and C in Figure 5. This is seen with more detail in Figure 6e, where the neighboring regions A and C are determined by the separatrices F 1 and F 3 , and in Figure 6f for region B, delimited by the separatrices F 6 and F 7 . Although regions A and C are close to each other, the double collision process is reversed when the impact parameter passes from the domain of one of them to the other. Thus, we find that while in regions A and B the atom first collides with the adsorbate and then with the flat surface, in the case of region C it is the opposite. The length of this latter region is nearly the same as the one for the shadow region S, given by Equation (23), because of the symmetry between the trajectories F 1 and F 7 On the other hand, regions A and B are defined, respectively, by the angular sectors ∆θ A = π/4 − θ i /2 − (sin) −1 (z 0 /a) and ∆θ B = π/4 − θ i /2, which are both dependent on the incidence angle.
Regarding trapping, although this simple potential function cannot lead to this phenomenon, it is interesting to note that, to some extent, F 3 and F 7 can be considered as permanently trapped trajectories, since they will keep evolving parallel to the surface (θ d = ±π/2). This is, of course, a rather weak case associated with a purely repulsive model, but still it is useful to understand in simple terms the types of dynamics that can be expected from a more refined classical model, such as the one specified by (2), concerning the presence of homologous pairs of trajectories (associated with single and double collisions) as well as the appearance of trapping.

Newtonian Level
The hard-wall model constitutes a sort of crude approach to the system studied, appropriate to explain some general features. However, a more realistic or refined model (closer to the experimental system) is going to include additional features, such as the rainbow phenomenon, which have also received much attention both experimentally [24,26] and theoretically [27][28][29], or the defect-mediated diffraction resonance [27,36], a kind of trapping induced by the presence of adsorbates on surfaces. These phenomena are associated to the particularities displayed by the potential model, such as the depth of potential wells, the stiffness of the repulsive region or the range of its attractive region. All these features are controlled by means of parameters that are found from best fit to the cross sections (intensities) experimentally measured; the optimal values eventually provide us with valuable information on the system analyze. In the particular case of the system discussed in this work, such information has to do with the way how the CO is attached to Pt(111) surface, which can be later used to better understand reactivity and diffusion properties.
The dynamics under the influence of the potential model (2) are illustrated in Figure 7a by means of a set of trajectories with initial conditions uniformly covering a wide range of impact parameters. Specifically, in this simulations b ∈ [−10.6, 10.6] Å, for E i = 10 meV and θ i = 20 • (again, as in the previous section, incidence out of the diagonal has also been chosen for these trajectories to stress some particular aspects). The attractive long-range term from Lennard-Jones contribution plays an important role here, since it is responsible for the permanent trapping (in this model) of He atoms along the Pt surface. This term accounts for the van der Waals interaction mediating between the neutral He atoms and the CO adsorbate, producing an effective transfer of energy from the perpendicular direction to the parallel one, such that the energy along this latter direction becomes larger than the incident energy to the expense of making negative the energy along the z direction. Of course, surface trapping is not totally permanent, since the presence of other adsorbates (not considered here, where we are working under single-adsorbate conditions) produces the opposite effect, that is, a trapped atom, after colliding with another neighboring adsorbates, may acquire enough energy along the z-direction (loosing it along the x direction) to escape from the surface.  (2) for an incidence energy E i = 10 meV. As in Figure 5, for a better illustration of all possible cases, also an incidence angle θ i = 20 • has been chosen. Moreover, the incident part of the trajectories has been represented with dashed line; (b) classical deflection function. While surface trapping gives rise to two kind of discontinuous regions, rainbow features manifest as local maxima (1) and minima (2); (c) Asymptotic-energy diagram. Here, trapping is detected through the two regions of the curve below the threshold E z = 0 meV, while rainbows manifest with two local minima (green squares; for rainbow 2, see enlargement of region A in the inset). Orange circles denote conditions leading to perpendicular deflection with respect to the flat surface.
A useful tool to systematize and analyze the different dynamical behaviors exhibited by the trajectories of Figure 7a is the deflection function, i.e., the representation of the deflection angle with which the He atoms asymptotically leave the substrate as a function of their impact parameter. This function is represented in Figure 7b and clearly shows that it is characterized by two types of regions. One of them is smooth and continuous, which means that the deflection angle increases or decreases gradually. The local maxima and minima within this type of region denote the presence of rainbows, i.e., deflection directions characterized by an extremely high intensity (leaving aside mirror reflection from the flat surface). In the figure, we observe the presence of two of these rainbows (denoted with the numbers 1 and 2), one of them nearly along the normal to the flat Pt surface. The other type of region is seemingly random, which is a signature of trapping-this makes the function in these regions to be time-dependent, since deflection will depend on the time at which it is computed. Nonetheless, this behavior is not to be misinterpreted with presence of chaos that we observe in analogous representations for He diffraction from corrugated surfaces [42]. The difference is that, while, in those cases, this random-like region has a fractal structure [42]; here, it is very regular. This can easily be seen by plotting more values of the impact parameter within this region. Then, as this region of the deflection function becomes more dense, we will be able to better appreciate its regularity and, therefore, the lack of an underlying chaotic dynamics.
To complement as well as to disambiguate the information provided the by deflection function, particularly in the trapping regions, an alternative representation can be obtained if we focus far from the adsorbate, where the system energy is separable, as mentioned above. Accordingly, unless there is an extra energy exchange because of the presence of other defects, the He atom will keep constant its energies along the x and the z directions (there is no coupling term between both degrees of freedom in the potential describing the flat Pt surface). This thus allows to consider energy diagrams, where the energy left along the z degree of freedom is compared to the total available energy, which coincides with the incident energy (E i ), and to the dissociation threshold (0 meV), which determines the minimum energy for unbound motion. This plot is shown in Figure 7c for E i = 10 meV. In this representation, we first note that a large portion of impact parameters gives rise to unbound motion, i.e., all the trajectories leave the Pt surface. Among these trajectories, those leaving the surface at an angle equal to a rainbow angle produce a local minimum in the energy diagram (green squares)-the absolute value of the rainbow angle can then be readily determined by means of the simple relation cos θ R = E R z /E i -while maxima (orange circles) correspond to trajectories that leave the surface perpendicularly (no energy along the x degree of freedom, although initially they all started with some energy in it). On the other hand, we also find some impact parameter regions for which the energy is negative. These regions correspond to conditions leading to trapping (oscillatory bound motion parallel to the flat surface), which are in correspondence with the discontinuous regions observed in the deflection function. Actually, from the energy diagram, we can determine with high accuracy the precise impact parameters that, at a certain incidence, will give determine the ranges of trapping. From the figure, we see that these limits are given by the intersection points of the energy curve with the zero-energy condition.
As we have seen, the energy diagram and the deflection function, in general terms, provide the same kind of information. Now, the energy diagram can be smartly used to determine pairs of homologous Newtonian trajectories as follows. Consider a given value for E z . All the impact parameters that are obtained from the intersection of a horizontal line at the selected value of E z with the energy diagram will provide sets of pairs of homologous trajectories. Some illustrative pairs are shown in Figure 8. In panel (a), for a selected energy such that E 1 > E z > 0, where E 1 is the energy (along the z direction) for rainbow 1, we find two pairs of trajectories, one back-scattered and another forward-scattered. Although distinguishing between single and double collisions is not as simple as with the hard-wall model, we still can perceive that within these two pairs of homologous trajectories one of them undergoes a single collision (denoted with black color), while the other shows something that could be related with a double collision (with blue color). If E z is below zero, we can still have two pairs of homologous trajectories, as seen in panel (b), although these trajectories exhibit permanent trapping. For energies above E 1 we can observe up to three or (for E z > E 2 ) four pairs of homologous trajectories, with analogous behaviors. Actually, unlike the hard-wall model, here we notice that there can be several pairs contributing to reflection symmetry interference, as happens in panels (c) and (d).
Finally, Figure 9 offers a comparison between the behavior displayed by a purely repulsive adsorbate (closer to the hard-wall model) and the attractive adsorbate in terms of the corresponding deflection functions for E i = 10 meV and normal incidence. The main difference between both models relies on a removal of the attractive term of the Lennard-Jones model. Notice that the repulsive model lacks the two rainbow features for ∆K R = ±1.95 Å −1 (θ R = ±26.50 • ), although both models keep nearly the same trapping rates, which has to do with the effective energy transfer from the perpendicular to the parallel directions due to the adsorbate curvature, and not that much with the presence of attractive basins around the adsorbate itself. Although not shown here, the same holds if the energy is increased. For instance, for E i = 40 meV and normal incidence, rainbows are observed for ∆K R = ±1.21 Å −1 (θ R = ±7.96 • ) with the full potential model (2), but there is a complete absence of them when the repulsive model is used instead.
where E 1 and E 2 are the energies (along the z direction) corresponding to rainbows 1 and 2 (see Figure 7b,c), respectively, and E i cos 2 θ i is the incidence energy (also along the z direction), with E i = 10 meV. Forward/backward deflected trajectories are denoted with solid/dashed line. Trajectories undergoing single/double collision/s are denoted with black/blue colors. All trajectories are started from a distance z i = 10.27 Å above the flat Pt surface (beyond z ≈ 6.35 Å, the interaction potential model (2) is negligible; see Figure 1a).  (4). Here, for simplicity, the incidence conditions are θ i = 0 • and E i = 10 meV. To compare with, the deflection function corresponding to the full interaction potential model (2) is also represented (blue solid line), obtained for the same incidence conditions.

Bohmian Level
In the literature there has always been a controversy concerning how the different diffraction features observed in intensity patterns should be assigned in the case of the atom diffraction by impurities on surfaces. In 1988, for instance, Yinnon et al. [27] suggested the use of the quantum flux, J, as an interpretational tool to understand trapping processes in this type of systems. A vector representation of this field provides a reliable representation of how the Fraunhofer, rainbow and trapping features arise, although is not unambiguous at all. This procedure was employed earlier on in reactive scattering by Wyatt [52][53][54]. Although some dynamical information can be extracted about the distribution and change of the flux, it is difficult to understand how the system eventually evolves. The last step in our journey is precisely the use of Bohmian trajectories to get an idea of what is going on this type of systems from at a fully quantum level and, beyond quantum flux based analyses, to understand the dynamics by causally connecting a point on the initial state with another point of the final state, following a well-defined trajectory in real time. This trajectory is obtained by integration of the equation of motion [12,14] often regarded as guidance equation. In this equation, formerly introduced by Bohm as a postulate [1,4], ρ and J are respectively the probability density and the quantum flux [55]; the velocity field v =ṙ is just the way how the probability density spreads through the configuration space in the form of a flux or current density. Regarding S, it is the phase field, which describes the local variations of the quantum phase and is typically obtained from the polar transformation of the wave function, The Bohmian trajectories shown below are obtained taking into account the value of the wave function at a given time (obtained by means of the propagator mentioned in Section 2), according to the second expression for the velocity field in Equation (24).
The analysis in this section is performed by studying the behavior displayed by sets of Bohmian trajectories with initial positions taken at a series of distances from the surface and uniformly distributed along the parallel direction for each one of those distances. Unlike the procedure followed in preceding sections, now the loss of translational symmetry caused by the presence of the isolated adsorbate on the Pt surface does not allow for studying the dynamics considering impact parameters only along the parallel direction for a given z value. This is the reason why different values of z are considered. This enables a better way to assign the different parts of the incident wave with final outgoing intensity peaks without ambiguity. These sets are displayed in Figures 10-12, with each set being labeled with the corresponding initial condition along the z direction referred to the center of the wave packet z 0 = 10.27 Å. According to these figures, quantum-mechanically the dynamics cannot be understood in the same local terms as in classical mechanics, where trajectories did not display a different behavior depending on their starting distance along the vertical direction with respect to the adsorbate. In order to obtain a complete knowledge of the diffraction process, trajectories have to be chosen from across the whole region covered by the initial wave function. Different regions will give rise to different diffraction features, but also different sets of trajectories will be able to probe the surface at a different level, being able to approach it very closely or, on the contrary, will bounce backwards far away, without even having touched it physically. In this regard, one wonders whether, contrary to what is commonly stated within the "Bohmian community", the (Bohmian) trajectories can be considered to reveal the "true" motion followed by the particles they are associated with-in the present case, the motion of individual He atoms. This is a rather challenging question (as well as metaphysical), which, to the author's best knowledge, has not still been unambiguously answered, that is, with a solid, irrefutable experimental proof. Therefore, taking on a pragmatic view, here Bohmian trajectories are considered as hydrodynamic streamlines that allow us to investigate the flow dynamics of the probability flux in configuration space and, therefore, to understand towards which directions atoms are more likely redirected (deflected) after being scattered from the substrate (without providing any particular information on how each individual atom really moves). Actually, to some extent, this random view of the atomic motion is the idea behind the work that Bohm developed in 1954 in collaboration with Vigier [56] and, later on, in 1989 with Hiley [57], or the approach developed in 1966 by Nelson [58] (although this latter approach has nothing to do with Bohmian mechanics, it introduces analogous stochastic concepts). In these examples, the quantum particle follows a random-like path as a consequence of the action of a sub-quantum random medium; when motions are averaged, particle statistics reproduce the results described by Schrödinger's equation and Bohmian trajectories arise as the averaged flow-lines associated with the solutions to this equation. Figure 10. Set of Bohmian trajectories with initial positions uniformly distributed along the x direction and fixed value along the z direction: z i = z 0 = 10.27 Å, which corresponds to the center of the incident wave packet (t = 0) above the clean Pt surface (beyond z ≈ 6.35 Å, the interaction potential model (2) is negligible; see Figure 1a). On the right-hand side, enlargement of the left plot near the adsorbate to illustrate the dynamical behavior displayed by the trajectories close to the substrate surface. Figure 11. Set of Bohmian trajectories with initial positions uniformly distributed along the x direction and fixed value along the z direction: (a) z 0 = z 0 − 3.18 Å; (b) z 0 = z 0 − 2.12 Å and (c) z 0 = z 0 − 1.06 Å, where z 0 = 10.27 Å, which corresponds to the center of the incident wave packet (t = 0) above the clean Pt surface (beyond z ≈ 6.35 Å, the interaction potential model (2) is negligible; see Figure 1a). In the corresponding lower panels, enlargement of the upper plots near the adsorbate to illustrate the dynamical behavior displayed by the trajectories close to the substrate surface. In Figures 11 and 12, sets of trajectories taken from below and from above the z-value selected in Figure 10 are shown. In the three cases displayed in Figure 11, we notice that those trajectories that start closer to the surface are going to contribute more importantly to the marginal portions of the outgoing wave (B i and C, as defined in Figure 3d). In sharp contrast, trajectories started in regions far from the adsorbate (see Figure 12) will contribute to the peaks related to small ∆K (A i ). This can be understood establishing a nice analogy with classical fluid dynamics and then introducing the notion of quantum pressure introduced by Takabayashi [59]. Accordingly, when the wave function is on the surface, we find that, while its lowest part is already bouncing backwards, the upper one is still moving downwards (towards the surface). In this situation, in the analogy of the wave function being associated with an ideal non-viscid and incompressible fluid, its upper part would be pushing the lowest one and then generating a remarkable pressure on it (this effect would increase with the incidence energy, although its duration would be shorter). This is something that cannot be seen directly on a wave-function representation, but that has an interesting counterpart in the case of the trajectories because we can see that those associated with the lowest part of the incident wave function are then pushed (or make evident the push) against the surface, remaining there for a rather long time. On the contrary, as the initial conditions are taken from upper regions of the initial wave function, the pressure on them will be lower and hence will bounce backwards further away from the physical surface. Thus, the evolution of these trajectories resembles in a closer way the behavior of classical trajectories, while the trajectories below are somehow forced to propagate parallel to the surface and escape by the borders of the wave.
The turbulent dynamics manifested by the trajectories close to the surface has to do with the web of nodal lines characterizing the wave function when it is on the surface. In particular, it is interesting that how sometimes the trajectories get trapped and whirl around some of these nodes, undergoing a sort of transient vortical trapping [32,33], different from the temporary trapping of the trajectories that will remain confined within the surface attractive well far from the adsorbate. The whirlpool motion displayed by such trajectories has actually an interesting property, namely the associated action is quantized, being an integer multiple (equal to the number of full loops) of a certain value. This type of motion is a confirmation of the former results found by Yinnon et al. [27]. Actually, recently, Efthymiopoulos et al. [60] formally determined the conditions under which such temporary nodal regions in scattering problems will appear, which happens to be in the regions where the amplitudes (modulus) of the incoming and outgoing waves become equal. This boundary, where both values are the same and that have important dynamical consequences, is what they called "separator".
Regarding the reflection symmetry interference phenomenon, this time it is not as simple as in the two previous cases because the trajectories obey the guiding rules imposed by the wave function and not the direct action of the potential. In the case of the rainbow, however, there is a more evident manifestation. Although Bohmian trajectories cannot cross one another at the same time, collectively they show a sort of rainbow-like precession, which is more prominent in the case of the trajectories shown in Figure 11: after the collision with the adsorbate, they start showing larger and larger deflections until getting trapped, then the start precessing backwards until reaching a maximum angle, and finally precess backwards again to smaller angles. In the case of the trajectories displayed in Figure 12, the behavior is analogous, even if in this case there is no trapping and the effect is more subtle. Thus, all the ensembles show a similar behavior, supporting the idea of the rainbow as being a sort of "global" effect, which translates into the "wings" observed in the intensity pattern shown in Figure 4.

Conclusions
In the last several years, Bohmian mechanics has been gaining ground as an appealing tool to deal with very different problems out of the area of the quantum foundations, its traditional environment. In such cases, the interest and relevance of Bohmian mechanics is emphasized by directly tackling a given problem with it. In the current work, however, the motivation has been a bit different. There are different trajectory-based approaches that have been or can be used to describe, analyze, understand and explain quantum systems and phenomena, even if they have different degrees of accuracy. The purpose here has been to establish an appropriate context to better understand the role of Bohmian trajectories within all those formulations as well as the kind of information conveyed by each one at its level of accuracy. To this end, we have considered a problem of interest out of the field of the quantum foundations, specifically the diffraction of helium atoms by a nearly flat platinum surface on top of which there is a carbon monoxide adsorbate. This is a problem that has been considered in the literature due to its intrinsic practical applications, although here it has been chosen due to its suitability to the purpose, since it is sufficiently simple to allow us its treatment at different levels. Accordingly, the system has first been analyzed with a usual wave-function-based framework, investigating the effects associated with two different He-CO/Pt(111) interaction potential models: • Hard-wall model. This model is in the form of an impenetrable (fully repulsive) wall, where the interaction is reduced to a sudden impact on the He atoms on the such a wall. The first model allows an exact asymptotic analytical treatment, convenient to elucidate the main mechanism observed in the diffraction pattern produced by single adsorbed particles on nearly flat surfaces, namely reflection symmetry interference. • Potential energy function. This interaction model is determined from fitting to the experimental data and constitutes a refinement of the previous one in the sense that there is detailed information on the intensity of the interaction between the incoming atom and the substrate at each point (in this regard, the hard-wall model is just a crude approximation). Thus, in spite of its lack of analyticity, unlike the hard-wall model, it provides us with a more realistic description of the diffraction process in real time, rendering information on additional physics, such as rainbow features or surface trapping.
Although at a different degree of accuracy, these two models provide us with explanation of the features observed in the experimental diffraction patterns. The question is how to interpret or understand these diffraction features or, in other words, to elucidate the physical mechanism responsible for each of such features (reflection symmetry interference, rainbows, or surface trapping). Typically, this is done by setting protocols based on quantum-classical correspondence, e.g., analyzing the system by means of classical trajectories and then comparing the results rendered by both the classical (trajectory) model and the quantum-mechanical (wave function) one.
Bearing that in mind, as has been seen in the preceding sections, here we have tackled the issue at three different levels: • Fermatian level. This first level is the simplest one, based on computing what has been here denoted as Fermatian trajectories, which are just the direct analog to optical rays reflected on a hard wall in a medium with constant refractive index. According to this trajectory model: -These trajectories have revealed that there are pairs of homologous trajectories, such that one of the peers undergoes single scattering off the interaction potential, while the other undergoes double scattering. The fact that a trajectory collides with the CO/Pt system at one point (single collision) or at two different points (double collision) is a function of the impact parameter. Accordingly, a simple mapping can be establish, which helps to easily localize regions of impact parameters that are going to produce homologous pairs of trajectories.

-
The mechanism of reflection symmetry interference is associated with these paired trajectories, which is explained in the same way that we explain interference from two coherent sources: interference maxima and minima arise depending on whether the path difference between the two paths (or virtual rays) joining each source with a given observation point on a distant screen is equal to an integer number of wavelengths or to half an integer, respectively. Although these paths are nonphysical (they are just a mathematical construct), they allow us to understand in simple terms the appearance of the alternating structure of bright and dark interference fringes. In the present case, the path length arises from the extra path length of the trajectory affected by the double collision with respect to the homologous pair with single collision.
-In addition, it has also been seen that two specific trajectories are deflected parallel to the surface, which can be interpreted as a mechanism precursor of the surface trapping mechanism that appears in more refined models, such as the Newtonian and the Bohmian ones.
• Newtonian level. On the next level, the Newtonian one, classical trajectories are obtained for the realistic potential energy surface describing the interaction between the He atoms and the substrate. In this case, it is not so simple to distinguish between single and double collisions, because the deflection of the trajectories near the surface, where the interaction between the He atoms and the CO/Pt surface is stronger, changes gradually very smoothly. However, we have been able to extract a series of interesting conclusions: -By means of an energy diagram (asymptotic energy along the z direction as a function of the impact parameter), we been able to devise a method that allows to determine in a simple fashion pairs of homologous (Newtonian) trajectories. This diagram is thus a suitable method to determine a behavioral mapping of initial conditions (impact parameters) for a given incidence direction (incident energy).
-Accordingly, also at this level, it is possible to find an underlying mechanism responsible for the reflection symmetry interference found in the corresponding quantum intensity patterns. Actually, interference patterns could be reconstructed in the same way as with the Fermatian model, although in this case we would be dealing with a space-dependent refractive index (the potential function) and the Newtonian trajectories would play the role of Feynman's paths. Nonetheless, although such a reconstruction is possible and the techniques are well known, this does not mean that trajectories, Fermatian or Newtonian, contain any information on the interference process; in both cases, they are only a tool to determine the interference pattern.
-Regarding the trapping phenomenon, it has been found to be more prominent, with an important amount of trajectories remaining trapped permanently along the surface. This is, however, only a temporary feature, since it may disappear as son as the trapped atoms find another adsorbate. In such a case, the collisions with this adsorbate may provoke an effective transfer of energy from the parallel to the normal direction, such that the will be able to eventually leave the surface.
-Finally, due to the attractive well surrounding the adsorbate, we have also observed the appearance of rainbow features, i.e., high accumulations of trajectories along particular deflection directions. However, rather than contributing with a specific, localized feature in the corresponding quantum intensity pattern, rainbows seem to manifest affecting them globally, i.e., giving rise to features that appear at different places. This has been noticed by computing exactly the same with an alternative repulsive adsorbate model, which lacks the surrounding attractive well and therefore does not give rise to the formation of rainbows.
• Bohmian level. The upper level here considered is the Bohmian one, where things change substantially if we note that the transition from the Fermatian level to the Newtonian one can be seen as a refinement associated with having a more accurate description of the interaction potential model, changing a hard wall by a "soft" wall. These are the main findings at this level: -First of all, since Bohmian trajectories are associated with a particular wave function, there is no freedom to choose a given set of initial conditions because depending on the positions selected relative to the region covered by the initial wave function, the trajectories are going to exhibit a different behavior. Thus, we have seen that while some of them are deflected quite far from the physical surface (more intense interaction region), other trajectories move just on top of it, displaying signatures of vorticality.
-To better understand that point, notice that Fermatian trajectories are only ruled by the law of reflection, while Newtonian trajectories are ruled by correlations between the two degrees of freedom, x and y, that can be locally established within the interaction region (i.e., the region where the interaction potential is stronger, near the substrate). In the case of Bohmian trajectories, the dynamics is not directly ruled by the interaction potential, but by a wave field that is able to (non-classically) convey information from everywhere in the configuration space (through its phase). This makes a substantial difference between classical (Fermatian or Newtonian) and Bohmian trajectories, which may lead us to think that direct comparisons or analogies must be taken with care. That is, nothing of what has been seen at the previous levels remains at the upper one, since it is not possible to form pairs of homologous trajectories.
-In this case, and contrary to the two previous models, the trajectories contain information about the interference process and, therefore, can be used to determine the fringe structure of the pattern by simply making statistics over them. If they are properly distributed across the region of the configuration space covered by the initial probability density, they will eventually distribute according to the final probability distribution by virtue of the continuity equation that they satisfy.
-Regarding rainbow features, present in the Newtonian model and also, with a weak precursor, in the Fermatian one, the only a similar behavior is observed, although it is difficult to establish a unique correspondence with the phenomenon of the two previous models. In the Bohmian case, taken the trajectories that start with the same value z 0 , it is seen that their final positions show, for some range of x 0 values, a certain "precession" as x 0 increases. However, it has not been possible to uniquely identify this phenomenon with the classical rainbow.
In the case of surface trapping, on the contrary, there same effect has been observed in the three models (again, in the Fermatian model it is only a weak precursor).
-Finally, it has also been observed that, depending on how close or far a Bohmian trajectory is started from the physical substrate surface, it will be able to reach this surface or just bounce backwards quite far from it (from what we could call an effective nonphysical surface). Actually, if the trajectories start close to the surface, they are influenced by the web of maxima developed (and sustained for some time) around the adsorbate, displaying a rich vortical dynamics.
To conclude, we can say that, although there is no one-to-one correspondence between classical and Bohmian trajectories, it is still possible to understand these two alternative descriptions in a complementary way, with one being the skeleton upon which the other rests, at least at a formal level. Classical trajectories have been and are used to understand in relatively simple terms why quantum distributions are as they are, in a way analogous to how an optical path allows us to understand and explain the appearance of wave phenomena, such as diffraction or interference. Bohmian trajectories are synthesized from wave amplitudes (wave functions), so the same underlying scheme should also be valid (i.e., using classical trajectories to understand why Bohmian trajectories evolve in the way they do). At the same time, Bohmian trajectories offer a clear picture of the evolution of the quantum system by monitoring the local evolution of the quantum flux, which provides some clues on dynamical aspects that otherwise would remain hidden (e.g., the development and effects of vortical dynamics, or the appearance of effective barriers).