{'fraction_non_alphanumeric': 0.03514522443776746, 'fraction_numerical': 0.012018574160065555, 'mean_word_length': 4.635454079014879, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-dimensional simulation of quantum reflection

A propagation method for the scattering of a quantum wave packet from a potential surface is presented. It is used to model the quantum reflection of single atoms from a corrugated (metallic) surface. Our numerical procedure works well in two spatial dimensions requiring only reasonable amounts of memory and computing time. The effects of the surface corrugation on the reflectivity are investigated via simulations with a paradigm potential. These indicate that our approach should allow for future tests of realistic, effective potentials obtained from theory in a quantitative comparison to experimental data.

While the effect on itself has been widely known since the early days of quantum mechanics [20], several papers, both theoretical [21][22][23][24] and experimental [3,5], were recently published that develop an insight into this purely quantum effect, which is crucial for modern surface science. Since most analytical treatments of quantum reflection rely on approximations, exact numerical efforts are needed to verify their regime of applicability. In parallel, matter-wave experiments enable tests of Casimir-Polder interactions [19,[25][26][27][28][29][30][31][32], that require quantitative simulations to bridge the gap between theory and experiment. A recent study investigates the effect of a periodically driven surface in one dimension [33] using a phenomenological atom-surface Casimir-Polder interaction potential. In two dimensions, a time-independent approach studies matter-wave diffraction due to quantum reflection from a doped surface [34].
While time-independent approaches are computationally more efficient when the mere final reflection amplitude is needed, a time-dependent method enables the treatment of time-dependent Hamiltonians, and offers more insight into dynamical details of the scattering process. In this work we apply for the first time a twodimensional time-dependent numerical scheme, that is suitable for the study of nonseparable potentials. We use a toy model potential, based on a phenomenological generalization of the well known 1/r 3 van der Waals potential [33], to investigate the zero-order effects of a periodic corrugation on the reflectivity and connect to the onedimensional (1D) results obtained in [33] for reference. Our results represent a proofof-principle that quantum reflection from realistic Casimir-Polder potentials can now be investigated quantitatively in two-dimensions (2D) via an optimized time-dependent numerical propagation procedure.

Numerical Method
For 2D dynamics, in which there is significant coupling between the two dimensions, semiclassical WKB-type methods, that work well in 1D static problems, are known to fail [21,[35][36][37]. Numerical wave packet propagation schemes, that are based on splitoperator techniques to treat the higher spatial dimensions (e.g. [38,39]) are not very stable either in the case of strong non-separability.
Therefore, the numerical procedure we choose makes use of the implicit, normpreserving Crank-Nicholson scheme with backward and forward substitution [40]. Assuming that the Hamiltonian is constant over a small enough time interval ∆t, H(x,ŷ), the finite difference time-dependent Schrödinger equation in Cayley's form reads: whereĤ is represented as an (N x ·N y )×(N x ·N y ) matrix. Equation (1) needs to be solved for ψ t+δt at each time step, so that the wave function can be known at any time during the scattering process. If three-point approximations for second derivatives are used along the x-and y-axes, the resulting structure of the Hamiltonian in the position basis is a tridiagonal block matrix with fringes [40]: The elements of the ψ(x, y) vector are indexed ψ(x, y) = [ψ 0,0 , ψ 0,1 , ...ψ 0,Ny , ψ 1,0 , ..., ψ Nx−1,Ny−1 ], where the subscripts represent the indices on the grid along the x-and y-axes respectively, and N x and N y are total numbers of grid points. It can be shown that this basis minimizes the matrix band to the width N y [41]. The Cayley form used in Eq. (1) has the advantage that a relatively large time step can be chosen without compromising the norm of the evolved state. Using a grid step which depends on position, at least in our realization with a three-point discretization, would make the corresponding matricesĤ non-symmetric, and consequently, we would lose the mentioned advantage of perfect norm preservation.
In order to optimize the solution to Eq.(1), the matrix which multiplies ψ t+δt is decomposed into the product of a lower triangular factor and its transpose, via Cholesky (or LL t ) decomposition [40], so that only one triangular matrix must be stored. For time-independent Hamiltonians, the Cholesky factor, which depends on the potential parameters and the spatial and temporal step sizes, only needs to be computed and stored once, assuming that these are kept fixed. Equation (1) can thus be solved by standard forward and backward substitution [42].
The band preservation property of Cholesky decomposition [43] implies that the resulting triangular Cholesky factors inherit the banded structure of the Hamiltonian, so that only a limited number of matrix elements of order N x × N 2 y needs to be saved. This reduces the memory requirements dramatically. Therefore, the memory cost scales linearly with the number of grid points along one (x), and quadratically along the other axis (y). This means that this method is most efficient for problems in which the numerical grids are very different along the two dimensions: In other words, for long and thin grids. In this method, periodic boundary conditions are easily implemented. This suggests to study systems with a periodicity along the axis, along which grid points are more expensive. In fact, it turns out that, in this basis, the additional matrix elements necessary to introduce periodic boundary conditions along the y-axis lie at the end of each x-block, within N y elements from the main diagonal. Therefore, they do not increase the size of the band, so that memory usage is unchanged and additional runtime is negligible.
Since we are interested in a scattering problem involving a particle, we assume for it an initial Gaussian wave packet with standard deviations σ x along the axis normal to the 1D surface (x) and σ y along the in-plane (y) direction. Its average position ( x , y ) and average momentum ( p x , p y ). The wave packet initially approaches the surface with p x < 0. The momentum representationψ(p x , p y ) is computed via a 2D Fast Fourier Transform at different timesteps [40], and the total reflectivity is best computed [33] by integrating the latter over the entire p y axis and the positive region of the p x axis, so that accounts for all reflected components of the wave packet at time t. Since R depends on time, its final value is taken only after it has reached the stationary regime. This may be obscured when, in addition to quantum reflection sharply oscillating spurious reflections arise from ψ(x, y) hitting the x-boundary of our numerical grid. We eliminate these by suppressing the transmitted part of the wave function by multiplying it at each timestep by a smooth sigmoidal filter function, whose parameters are optimised to minimise such artificial reflections (see [33] for details).

Paradigm Potential
In one dimension, the atom-surface interaction can be described, in the short-range v/d Waals limit as: where C 3 and r 0 are material coefficients and r is the distance of the atom along the normal to the point-like surface [33,44,45]. In order to simulate the effect of a small corrugation along the surface, we generalise equation (3) by introducing a periodic modulation in the distance to the surface, which leads to a first order approximation to the results in [10,46], such that where A, L and φ are the amplitude, wavelength, and phase of the corrugation respectively. This defines our paradigm potential for testing the influence of the periodic corrugation on the quantum reflectivity R(t → ∞). The singularity, which occurs as x → A sin ( 2π L y + φ), is removed by introducing an artificial cutoff length ∆, beyond which the potential is continued as a parabola with vertex in r(x, y) = 0. For r(x, y) < 0, we impose ∂V ∂x = 0, so that V (x, y) becomes independent of x. Continuity and differentiability are imposed at the cutoff, as a generalisation of the approach introduced by [33] for the 1D case.
The complete, smooth 2D paradigm potential can thus be written as Introducing a cutoff length ∆ for the potential induces oscillations in R as ∆ is varied, as shown in later sections, and the effective reflectivityR can be computed by averaging over these oscillations [33]. We end this section by noting that this continuation is effectively 1D, as it simply shifts the point at which the parabola is matched with Eq. (5), so that the locus of points at which the potential is continued satisfies the equation x − A sin ( 2π L y + φ) = ∆, meaning that it is an equipotential line. Our simplified potential is invariant under a translation which leaves r(x, y) in Eq. (4) unchanged. This enables such a simple continuation to be straightforwardly implemented while preserving the 2D nature of the problem. In fact, an additional arbitrary parameter would need to be introduced in the continuation if less symmetric potentials are used. Physically, a continuation along the x axis is justified by the fact that the gradient of the potential, which is the key physical parameter of this problem, has a y-component which is proportional to the ratio A/L, so that in the small amplitude regime the steepest potential gradient is always along the x-axis.

Results of numerical simulations
We now present the results of our tests on convergence and further investigations of the effects of a small corrugation on the reflectivity for the potential defined in Eq. (6). We consider a corrugation wavelength L = 100 nm, which coincides with the size of our numerical grid along the y-axis, for which periodic boundary conditions are used. Orthogonal to it our grid spans the region −1.5 µm ≤ x ≤ 5.0 µm along the surface normal. The interaction constant was chosen C 3 = 4.0 × 10 −50 J which corresponds to the v/d Waals potential between a 3 He atom and a Au surface [44,45] (m ≈ 5.01 × 10 −27 kg [47]). A sketch of the potential with the initial wavepacket is shown as a 3D plot in figure 1.
The departure from the flat surface case is investigated by gradually increasing the amplitude A of the corrugation, which introduces a coupling between the two dimensions in the Hamiltonian. The width of the incoming Gaussian wave packet along the normal axis x is fixed at σ x = 80 nm. The initial wave packet is always placed at x(0) = 2.0 µm, with an initial average velocity ( v x (0) , v y (0) ) = ( p x (0) /m, p y (0) /m) = (2.0, 0.0) m s −1 .

Convergence
Timestep convergence tests confirm that dt ≤ 6.0 ns guarantees good convergence, with relative residuals of order 10 −4 for all spatial step sizes at which spatial convergence is achieved. Therefore, in the following simulations, we use dt = 5.0 ns. To test our numerical machinery, we study convergence with respect to spatial grid sampling. Simulations were performed using a small corrugation A/L = 0.1. We fix the phase φ of the corrugation to zero at the initial position of the wave packet y = 0 . The width of the initial wave packet along the y-axis is σ y = 8 nm. The results of these tests w.r.t. cut-off distance ∆ are shown in figure 2. In each subplot, the reflectivity is plotted as a function of the number of grid points used along the normal axis N x . Powers of two are used for the number of grid points in order to facilitate the use of a Fast-Fourier Transform [40]. The different symbols represent the level of precision along the y-axis, N y .
Each row in fig 1 corresponds to a fixed value of the cutoff length ∆. The left column of plots show R in logarithmic scale, whereas the right one shows a zoom into the convergence region on a linear scale. Figure 1 shows how smaller values of ∆ require a higher spatial resolution for convergence. This is expected for the toy potential of Eqs. (5) and (6): Since the potential becomes singular as x → A sin (ωy + φ), the large potential gradient encountered at small values of ∆ results in a poor sampling of the potential landscape. As a result, a higher precision in the spatial grid is needed for convergence. The abrupt increase of R as the numerics breaks down is a known effect from previous studies [33,41]. At that point, the potential is effectively sensed by the numerics as a step, which generally has a higher reflectivity than a soft potential [2]. The relative errors associated to the present numerical method, which was thoroughly tested in 1D and compared to exact solutions [33], are of the order of ∆R/R ≈ 1%. The aforementioned 1D treatments show that convergence with respect to ∆ is achieved when ∆ is in the nm regime and fractions thereof, where it can be extracted efficiently via logarithmic averaging [33]. Taking this into account, we observe that good convergence can be achieved with a number of grid points along the x-axis N x ≥ 2 15 and along the y-axis N y ≥ 2 7 up to a cut-off of ∆ = 3 nm, whereas a finer grid enables convergence at ∆ = 1.5 nm. Finer grids can still be used, especially along the x axis, in order to reduce ∆ to the sub-nanometre regime.
The scaling of the computational costs with different grid sizes is shown in tables 1a through 1c.

Effects of cutoff and corrugation
Next, we systematically study the effect of the cutoff length ∆ on the reflectivity. Starting from a A/L = 0 (i.e. the flat surface), we increase the corrugation amplitude gradually. Figure 3(a) shows the oscillations in R as a function of ∆, which appear analogous to those in our 1D results [33]. The reflectivity for the flat surface is compared to 1D time-independent simulations (continuous line) for different cutoff lengths ∆. N y = 2 7 was used in all simulations. For smaller values of the cutoff length, the symbols for N x = 2 14 and N x = 2 15 can be distinguished, which signifies poor convergence. The 1D results are accurately reproduced by our 2D method in the flat surface case; we find a relative difference between the two methods of less than ∆R/R ≤ 1% down to ∆ = 3.0 nm; below this cutoff convergence is lost. We also calculate the reflectivity as a function of A/L for fixed cutoff points (figure 3b). We observe that R changes with the corrugation amplitude A for all values of the cutoff. These deviations are smaller for larger values of the cutoff. This is expected, because if the potential is parabolically continued very far from the surface, the effects of the corrugation are too small to influence the dynamics significantly.

Effect of the width of the wave packet
In order to observe the effects of coupling introduced by the corrugation, we study the effect of increasing the width of the incoming wave packet σ y on the calculated reflectivity. Since we are only interested in the effect of widening the initial wave, we  The computational cost of the simulations for different numbers of grid points along the x (columns) and y (rows) axes. Table 1a shows the memory usage in GB, 1b the approximate number of hours needed to perform a single Cholesky factorisation (LL T ), and 1c gives an upper bound for the propagation time per single time step (dt=5 ns) in seconds. Machines used for computations requiring more than 4GB of memory feature AMD Opteron 6282 SE CPUs. just use the part of the Gaussian, with width σ y in the second (y) direction, which is within the considered interval L, with proper normalization of the y marginal of the modulus of the wave function. For σ y L, much smaller than the periodicity L in y direction, this corresponds to a usual Gaussian, whilst in the other limit σ y L we obtain an essentially flat profile in the y direction. The results are shown in figure 4 for N x = 2 14 and N y = 2 7 . For wider wave packets, an increase in the reflectivity is observed for the corrugated surface (A/L = 0.1), whereas essentially no dependence is observed for the flat surface. This is exactly what is expected: for the flat surface no scattering occurs in the y-direction, whereas the corrugation couples the x and y momentum components of the wave packet. Since the initial state is a Gaussian wave packet with v y = 0 m s −1 , a wider distribution in position space (along y) results in a momentum distribution more localised around p y = 0 kg m s −1 . Because quantum reflection decreases steeply with increasing incident momentum [33], waves having fewer spectral components of larger momentum are reflected more. As the width of the wave packet approaches the size of the numerical grid, the reflectivity saturates, so that the incoming wave packet effectively becomes uniform along the y-axis. For σ y < 4.0 nm, the grid sampling is not sufficient for a precise resolution of the initial wave packet, and a sharp, non-physical decrease in the reflectivity is obtained (for very small σ y ).

Conclusions and Perspectives
In this work, we have presented numerical simulations on quantum reflection, using an optimised Crank-Nicholson scheme in combination with a Cholesky decomposition tailored to a banded matrix system. We demonstrate that this 2D approach reproduces our earlier 1D results in the limit of zero corrugation. In addition, we show that this numerical method can be applied for the time-dependent simulation of quantum reflection from a realistic 2D non-separable potential. While this method is especially suitable for periodic problems such as uniaxially corrugated surfaces, it can in principle be used for non-periodic structures as well (modulo additional computational costs). As a first practical result we have shown that the coupling between the two dimensions introduces a dependence of the reflectivity on the corrugation amplitude.
While the feasibility of this type of computations has been demonstrated here for a phenomenological 2D non-separable potential, the possible parameter space is still very broad and can be explored and adjusted according to the specific needs of the experiments. We postpone such investigations to the future when more realistic 2D potentials will be available for modeling actual experiments. Our numerical technique allows us to investigate quantum reflection from any given conservative potential and can test theories that produce different effective potentials, amongst each other, and with respect to experimental results.
As shown in [33], higher order approximations for second derivatives in the Hamiltonian enable better precision. Since the consequent extra terms along the y-axis do not increase the size of the band, the computational costs for their implementation should be negligible, thus obviating for the extra cost of grid points along the y-axis (quadratic in N y ). In addition, the time step dt might be made adaptive by performing Cholesky factorisation at different points during the propagation, depending on the numerical precision needs.
From a numerical point of view, the use of the Demko-Moss-Smith theorem for the decay of matrix elements in the inverse of a banded matrix away from the band [48], may save memory since only a smaller approximation of the inverse of the matrix in the lhs. of Eq. (1) is needed. This may allow us to replace the iterative forward/backward substitution procedure by matrix-vector products only, which are simple to implement and fast to execute. Alternatively, spectral approaches may be considered, which exploit the symmetry and structure of the Hamiltonian, such as the one proposed in [49]. Finally, the forward-backward sutstitution procedure might be implemented in parallel, thus reducing dramatically the time required for each propagation step.
We end this paper by discussing some practical consequences of our work: the length scale for quantum reflection ideally matches that of the periodicity and corrugation height encountered in nano-materials. For conventional atom-surface scattering longrange order and atomic flatness are required not to loose the entire beam intensity in diffuse scattering. In quantum reflection, however, the turning points lie much further away from the surface, which makes our numerical propagation method well suited even for large roughness that nano-materials typically exhibit. Here, not only their static properties, but also their dynamics can be investigated, i.e. when the Hamiltonian is explicitly time-dependent. Moreover, our method allows for studies of time-dependent details of the quantum-reflection scattering process itself that cannot be addressed directly using time-independent approaches. As an example, the propagation of surface plasmons and bulk plasmonic excitations that extend all the way to the surface are very relevant processes in organic electronic materials. These open questions can now be investigated systematically. Furthermore the dynamics within the novel magnetic materials developed in recent years can studied with quantum reflection, using scatting particles with a spin (such as 3 He used in atomic beam spin echo spectrometers). An example of the search for structure and dynamics in micro-structured artificial magnetic assemblies is discussed in [50]. Note that these magnetic structures do not necessarily possess a physical corrugation that comes along with the magnetic one. The ever decreasing physical size of a single bit in magnetic storage devices and, therewith correlated, their increasing volatility, offer a true playground for the time and length scales of our numerical method. From a more fundamental point of view, our method may steer and optimize the design of a (dynamic) beam splitter for atoms that functions without laser light.