The discrete Green's function method for wave packet expansion via the free Schr\"odinger equation

We consider the expansion of wave packets governed by the free Schr\"odinger equation. This seemingly simple task plays an important role in simulations of various quantum experiments and in particular in the field of matter-wave interferometry. The initial tight confinement results in a very fast expansion of the wave function at later times which significantly complicates an efficient and precise numerical evaluation. In many practical cases the expansion time is too short for the validity of the stationary phase approximation and too long for an efficient application of Fourier collocation-based methods. We develop an alternative method based on a discretization of the free-particle propagator. This simple approach yields highly accurate results which readily follows from the exceptionally fast convergence of the trapezoidal rule approximation of integrals involving smooth and rapidly decaying functions. We discuss and analyze our approach in detail and demonstrate how to estimate the numerical error in the one-dimensional setting. Furthermore, we show that by exploiting the separability of the Green's function, the numerical effort of the multi-dimensional approximation is considerably reduced. Our method is very fast, highly accurate, and easy to implement on modern hardware.


Introduction
The free expansion of wave functions describing massive, non-relativistic particles is an important computational problem because of its high relevance for a number of applications in matterwave interferometry [1,2], where atoms are cooled down to temperatures in the microkelvin or even nanokelvin range, so that their motion becomes essentially quantum and thus dominated by wave phenomena, such as interference.
Initially, the atoms are confined in a magnetic or optical trap undergoing coherent manipulations like beam-splitting and beam-recombination operations.To perform a measurement, the atoms are released from the trap and their interference pattern is detected after a certain time of flight in the field of gravity [3,4], which can be accounted for using an accelerated reference frame.
If the measurement is performed not with individual atoms or ions [5], but with Bose-Einstein condensates consisting of many thousands of atoms [6,7], atom-atom interactions are not negligible during the short initial stage of the expansion.However, this can be easily taken into account by solving the nonlinear Gross-Pitaevskii equation [8,9] for a couple of milliseconds on a moderately extended spatial grid until the largest part of the interaction energy has been converted into kinetic energy.The subsequent expansion is basically ballistic and the problem of the time-of-flight expansion boils down to solving the free Schrödinger equation in three dimensions.
This two-stage approach has been successfully applied in our recent three-dimensional simulations [10] of a bosonic Josephson-junction, where we were able to reproduce the experimental results in [11] very well.This problem looks simple, but, in fact, it is not.First of all, the expansion time is limited by the size of the laboratory setup and, therefore, not long enough to ensure the applicability of the stationary phase principle which yields an asymptotic expression for the wave function in the far-field limit 1 .On the other hand, due to the tight external confinement and the Heisenberg uncertainty principle, a considerable amount of energy is stored in the trapped quantum gas.This energy converts into kinetic energy when the atoms are released from the trap leading to a fast expansion of the atomic cloud in the tightly confined directions.Depending on the strength of the confinement and the time of flight, the volume of the atomic cloud increases by a very large factor (up to ∼ 10 3 ).
Previously the free expansion phase was handled by a simple Fourier collocation method consuming excessive amounts of computational resources.Imposing periodic boundary conditions, the free wave packet expansion problem can be solved using the time splitting spectral [12] (Fourier split step) method.To this end, the wave function in the free Schrödinger equation is replaced by a trigonometric polynomial and the equation is required to hold at the collocation points x j , j = 0, . . ., J. Since the potential in the free Schrödinger equation is zero, the method reduces to a single time step which is computed in O(J log J) time using the fast Fourier transform (FFT).The resulting numerical procedure is what we refer to as the Fourier collocation method.Unfortunately, due to the vast expansion of the wave function, the size of the computational domain and therefore the number of collocation points J is required to be very large 2 .In two and in particular in three spatial dimensions the number of required grid points becomes astronomically high.Eventually, the number of required grid points is so large that the numerical approximations of the initial and the expanded wave function cannot even be represented in local memory.Furthermore, despite the favorable complexity of the FFT-algorithm the computational effort in three spatial dimensions is considerable, revealing the need for a better numerical procedure.
The memory problem described above has been addressed in [13].The method uses two different spatial grids to represent the initial wave function ψ 0 on Ω 0 and the final expanded wave function ψ on Ω.Both grids employ J grid points but the grid spacing of the final grid is enlarged by a factor m ∈ N. Consequently, the memory requirements are reduced dramatically which, however, comes at the numerical costs of m applications of the FFT of size J.In fact, the algorithm in [13] can be seen as a clever way of computing only every mth value of the numerical approximation in the Fourier collocation method.
The idea of computing an approximation of the expanded wave function on a much coarser spatial grid is motivated by a simple observation.Wave functions are non-observable quantities.The detector measures essentially the density, i.e., the square of the absolute value of the wave function.If we are only interested in the density, there is no need to resolve the fine details in the real and imaginary part of the expanded wave function.Nonetheless, also the method presented in [13] is based on periodic boundary conditions and therefore the final domain Ω is still required to capture the entire non-zero part of the expanded wave function.In other words, the size of Ω is determined by the fastest moving parts in the initial wave packet.If Ω is too small, parts of the wave packet will contaminate the numerical solution by periodically reentering the domain from the boundaries.
Alternatively, one might consider the application of a domain truncation technique like complex absorbing potentials [14,15], perfectly matched layers [16,17,18,19] or the recently introduced Fourier contour deformation approach [20,21].However, this idea is not very helpful in solving the wave packet expansion problem since we are particularly interested in the interference pattern forming at late times which implies large spatial scales.In other words, any domain truncation technique would literally eliminate most or all valuable information accessible only in the expanded interference pattern.
An effective way to get rid of the above mentioned boundary condition issues is to consider the integral formulation of the solution.In particular, we propose to employ a simple discretization of the free particle propagator of the free Schrödinger equation.This approach to solve the free wave packet expansion problem seems so obvious that it is hard to believe that it has not been used before.One reason for this could be that a direct discretization of the single particle Green's function seems to be too simple to yield accurate results.Another reason might be that it was believed that the numerical effort to evaluate the discrete free particle propagator is quadratic in the number of grid points J even in spatial dimensions higher than one.It turns out that none of these assumptions are true.
In fact, we show that the most simple discretization of the underlying convolution formula yields stunningly accurate results.In all examples presented below only a very modest number of grid points J is needed until the numerical error hits the inevitable barrier caused by rounding errors in the double precision arithmetic3 .However, the spectacular convergence rate observed in the examples is a well-known effect in the numerical analysis literature and in particular in the field of pseudospectral methods.It is based on the fact that the trapezoidal rule4 approximation of an integral for a rapidly decaying and sufficiently smooth function converges at a high-order algebraic, spectral, exponential or even super-exponential rate with respect to the number of employed discretization points.In the context of this magical phenomenon we would like to mention the pioneering work in [22] as well as the famous review in [23].
With regard to the second issue, we suspect that the separability of the multi-dimensional problem has been overlooked.Quite obviously, the corresponding d-dimensional Green's function can be factorized into d one-dimensional free particle propagators.By exploiting this simple observation the numerical effort of the multi-dimensional discrete Green's function approximation is reduced tremendously.The numerical effort to solve the three-dimensional problem is, for example, no longer in O(J 2 ) but in fact only in O(J 4/3 ).Here, for convenience only, we have assumed that the number of grid points J = J 1 J 2 J 3 needed to discretize the initial wave function ψ 0 on Ω 0 coincides approximately with the number of grid points K = K 1 K 2 K 3 used to approximate the expanded wave function ψ on Ω.In general, this assumption is not needed.It is rather possible to employ two separate spatial grids wherein J 1 , . . ., J d and K 1 , . . ., K d are adapted to the problem at hand.Moreover, unlike in the case of the Fourier collocation method, the spatial grid corresponding to the final approximation is not required to cover the entire non-zero part of the expanded wave function.In fact, it is possible to consider any finite d-dimensional rectangular domain allowing to investigate the most interesting part of the wave function only.
The article is organized as follows.In the remaining part of this introduction we will introduce the main problem in mathematical terms and give an important one-dimensional example.By means of this example we also demonstrate the limitations of the widely used stationary phase approximation.In Section 2, we introduce the one-dimensional discrete Green's function approximation including an error analysis for two important classes of initial wave functions.Finally, Section 3 covers the discrete Green's function approximation for the multi-dimensional problem.In particular, we present an implementation of the approximation in three spatial dimensions which is then used to solve another set of non-trivial examples.

The free wave packet expansion problem
We consider the free Schrödinger equation with the boundary condition lim or rapidly decaying for |x| → ∞.We are interested in computing a numerical approximation of ψ = ψ(•, t) on at a fixed final time t > 0. From the physics point of view, Ω corresponds to the area accessible by the imaging system in an experiment.Typically, the initial wave function is expected to expand along all coordinate axes, and hence, Ω 0 ⊆ Ω is a reasonable requirement.However, we will see that this requirement is not needed and that the ability to compute ψ on a finite but otherwise arbitrary rectangular domain is a valuable feature.

Scaling
In the numerical experiments presented below we measure length in units of 0 = 1 × 10 −6 m, mass in units of the 87 Rb atom mass m 0 = 87 amu ≈ 1.45 × 10 −25 kg and time in units of t 0 = m 0 2 0 / ≈ 1.37×10 −3 s.As a result of this scaling we have = m = 1 which is employed throughout all equations presented below.We note that this scaling is also used in ultracold atom experiments and hence the numerical simulations are close to real-world observations.

Example: Superposition of two phase-shifted Gaussians
As an illustration we consider a simple but important example.The initial wave function is given by a superposition of two phase-shifted one-dimensional Gaussian wave packets where σ = 1/2 and δ = 5/2.We note that ψ 0 is not normalized to one.Due to the linearity of the free Schrödinger equation, the normalization is actually irrelevant.The corresponding exact solution to the free Schrödinger equation is given by [24] ψ(x, t) = 1 1 + i(t/τ ) (3) for x ∈ R and t ≥ 0 using τ = 2σ 2 .

Asymptotic time evolution
The solution to the one-dimensional free Schrödinger equation may be written as where denotes the Fourier transform of an integrable function f : R → C. Using Eq. (4) reads which by means of the stationary phase principle [25] yields the approximation The exact order of convergence is a little tricky to calculate as the phase factor S x,t (ξ) depends on the parameter t itself.However, replacing ψ 0 in (5) with the initial wave function (2) of the example in the previous section yields where we substitute ξ = x/t.The real and imaginary part of the asymptotic approximation (6) and the exact solution (3) are shown in Fig. 2  ) and is still relevant for expansion times 50 ms in a real experiment.

Discrete Green's function approximation for the one-dimensional problem
Figure 3: Real and imaginary part of the one-dimensional kernel G (1) in ( 8) for t = 1, t = 2 and t = 4.
Since the solution of the multi-dimensional problem can be reduced to the solution of several one-dimensional problems, we consider the one-dimensional problem first.

Discrete convolution
The Green's function formalism for the one-dimensional free Schrödinger equation reads where is the one-dimensional free-particle propagator [26].Fig. 2 shows the real and imaginary part of G (1) for three different final times t.
The initial wave function ψ 0 is assumed to decay rapidly.Alternatively, ψ 0 is assumed to be compactly supported on Ω 0 = [−L/2, L/2] for some L > 0. In either case we let denote a discrete representation of the initial wave function at the grid points for an even5 integer J ∈ N.
Our aim is to compute a numerical approximation of ψ(•, t) on the interval Ω = [a, b].To this end, we introduce the grid points and the approximation ψ = ψ0 , . . ., ψK−1 , where By means of ( 7) we obtain The integral is then replaced by the discrete convolution In practice, all approximations are computed simultaneously using with the discrete propagator Remark 1.In the further course of this paper we will frequently consider the error ψ − ψ ∞ and the relative error denotes the exact solution and ψ ∈ C K is the approximation in (12).We point out that the accuracy of each of the K approximations ψk , k = 0, ..., K − 1 is independent from the parameter K but depends only on the number of quadrature points J.The parameter K, on the other hand, determines the resolution of the expanded wave function.
In the examples below we use relatively large values of K which allow for a nice visualization of the computed approximations.

Application to the example from Section 1.3
As an example, we apply the discrete Green's function approximation (12) to the initial wave function in (2).In particular, we compute numerical approximations of the expanded wave function ψ(•, t) on the intervals Ω = [ −20, 20], Ω = [−40, 40] and Ω = [−80, 80] corresponding to the final times t = 2, t = 8 and t = 32, respectively.We always choose L = 20 for the length of the initial domain Ω 0 = [−L/2, L/2] and K = 1024 for the number of grid points in the final domain Ω.
Fig. 4 shows the relative error as a function of the number of grid points J for the three final times t = 2, t = 8 and t = 32.The results are presented for different values for the widths σ of the Gaussian wave packets in the initial wave function (2).From the shape of the curves and the fact that we employ a semi-logarithmic scaling it is clearly visible that the numerically observed convergence rate is faster than exponential.Only a very modest number of grid points J is needed until the relative error hits the inevitable barrier caused by rounding errors in the double precision arithmetic.
As expected, the wave packet expansion problem becomes more difficult to solve when the smallest feature size (here the width σ) in the initial wave function gets smaller.The same holds true if the final time t becomes smaller which is a direct consequence of the fact that the propagator in (8) becomes singular for t → 0. Extremely small final times t are, however, not of practical relevance in the free wave packet expansion problem7 .

Error analysis for analytic initial wave functions
In this section we consider rapidly decaying initial wave functions which are defined on the real line R and have an analytic extension to the strip in the complex plane for some c > 0.
Our aim is to estimate the error of the discrete Green's function approximation (12) at a fixed final time t > 0. To this end, we let with x = L/J and x j = (j − J/2) x denote the numerical approximation of the solution Within the Green's function approximation (12), the expression in ( 14) is evaluated at x = x k , k = 0, . . ., K − 1 using x k = a + k x and x = (b − a)/K.Since all grid points are located inside the interval [a, b], the error is bounded by In order to evaluate (15) we first consider the error and h = x we may write where we introduced the expressions and The first expression is the truncation error which in our application is independent from x. Since ψ 0 is assumed to be a rapidly decaying function, the truncation error in a typical application is extremely small.The second term is the discretization error The discretization error implicitly depends on x and t.It decreases exponentially fast, provided f meets the requirements of the following theorem [23]: Theorem 1.Let f be a complex function defined on the whole real line.Suppose further that f has an analytic extension to the strip Z c in (13) for some c > 0 and f (z ) → 0 uniformly as |z | → ∞ in the strip.Moreover, for some M , it satisfies for all y ∈ (−c, c).Then, for any h > 0, exists and satisfies Moreover, the quantity 2M in the numerator is as small as possible.
Provided the truncation error can be neglected, Theorem 1 yields some constant for some constants α, β > 0 and hence the error decrease exponentially fast (or faster) with the number of grid points J.
For illustrative purposes, we consider the initial wave function ψ 0 in (2) again.In particular, we demonstrate the calculation of error bounds for the approximation of the expanded wave function at the final time t = 8.We use the same set of parameters that was used to compute the approximation depicted in the third column of Fig. 1.Estimating the truncation error using (19) as well as calculating M * via Theorem 1 is a relatively simple but tedious exercise that is shown in Appendix A. The final result of these calculations is given as follows: We note that the initial wave function in the example is an entire function, meaning that it has an analytic extension to the whole complex plane C. Since G (1) in ( 8) for fixed t > 0 is an entire function as well, the integrand f in (17) has also an analytic extension to the whole complex plane C. In particular it has an analytic extension to the strip Z c in (13) using any positive number c.This is the reason why the error bound in (23) is in fact true for any c > 0.
In Fig. 5 we show the error bounds Γ c (J) for c = c using c = 0.2 and = 1, . . ., 30.We also show the maximum difference of the numerical approximation ψ ∈ R K and the exact solution ψ ∈ R K using K = 1024 grid points in the final domain Ω = [−40, 40].It is immediately apparent that the envelope of the calculated error bounds is only slightly larger than the errors of the actual approximations.

Error analysis for compactly supported initial wave functions
We now consider initial wave functions of the form where ) for some m ≥ 0. Like in the previous section our aim is to estimate the error of the discrete Green's function approximation (12) at a fixed final time t > 0. Using for fixed x ∈ [a, b].To simplify the notation, we introduce the function and set h = x .Moreover, we include the rightmost grid point x J = L/2 into the set of grid points (9).Since u(x 0 ) = u(−L/2) = 0 and u(x J ) = u(L/2) = 0 we also find f (x 0 ) = f (−L/2) = 0 and f (x J ) = f (L/2) = 0. Consequently, we may write which coincides exactly with the trapezoidal rule approximation of the integral Estimates for the error of the trapezoidal rule approximation are readily available.Based on the Euler-Maclaurin formula [27,28,29] they are typically formulated for real-valued functions.The case of a complex-valued function requires only a minor modification and the corresponding result is given as follows: for some m ≥ 0. Further define h = L/J, x j = −L/2 + jh, j = 0, 1, . . ., J for some J ≥ 1 and let denote the trapezoidal rule approximation of the integral Then, the error is bounded by where and The factors B denote the Bernoulli numbers for ∈ N 0 .
Using Theorem (2) we can estimate the error (25) for a fixed x ∈ [a, b].Moreover, it is possible to compute error bounds C m (J) such that for every integer J > 0. Since the calculations are of very technical nature we leave the details of our approach to the interested reader, see Appendix B.
As an example, we consider the compactly supported initial wave function ψ 0 in (24) using with the polynomials respectively.Densities, real and imaginary parts of these approximations as well as the initial condition are depicted in Fig. 6.The computation of an exact reference solution should, in principle, be possible but the final expression would look incredibly complicated.Nonetheless, by means of Theorem 2 we are able to demonstrate that the numerical approximations shown in Fig. 6 have been computed with utmost precision.
The fact that u is a polynomial with vanishing boundary values implies f ∈ C ∞ 0 ([−L/2, L/2]).We are therefore free to choose any m ∈ N 0 in Theorem 2. For t = 8 the obtained error bounds C m (J) are shown in Fig. 7 using J = 16, . . ., 512 and m = 0, . . ., 10. Quite obviously, the guaranteed speed of convergence is not as fast as in the previous examples.It should however be noted that our estimates are comparatively coarse as they are based on a simple expansion of the derivatives of the integrand.Nonetheless, using J ≈ 260 grid points the error is guaranteed to be smaller than the machine precision which still represents a remarkable result.Furthermore, we find that C m (J) decrease like O(J −10 ) if m and J are sufficiently large.This convergence behavior can be explained by the fact that u (n) (±L/2) = 0 for n = 1, 3, 5, 7 but u (9) (±L/2) = 0.The same applies to the Figure 7: Error bounds Cm(J), m = 0, . . ., 10 in (28) for the discrete Green's function approximation in case of a compactly supported initial wave function ( 24) using the high-order polynomial u in (29).

Numerical effort of the one-dimensional approximation
The numerical effort involved in the evaluation and application of G in ( 12) is in O(JK) or O(J 2 ) if J and K are comparable in size.If Ω 0 = Ω, or more precisely, if J = K and the same operation requires only O(J log(J)) elementary numerical operations.In that case the discrete free particle propagator G is a square matrix where each descending diagonal of G from left to right is constant.Hence, G is a Toeplitz matrix and by embedding G in a circulant matrix of size 2J, the matrix-vector product in (12a) can be computed as follows [30]: First, we define the two vectors Finally, the last J components of w are discarded which yields where x = x = L/J.The complexity of the above method is in O(J log(J)) provided the DFTs are computed using the FFT algorithm.Realistically, however, the condition in (30) represents a major restriction.Like in the Fourier collocation method, the same set of grid points is used to approximate the initial as well as the expanded wave function and therefore the number of required grid points J in a typical wave packet expansion problem is very large.This is in strong contrast to the discrete Green's function approximation in its most simple form (12) which allows to compute an approximation of the expanded wave function on an arbitrary interval Ω = [a, b] ⊂ R using a customized number of grid points K. Taking into account this additional flexibility, the method in (12) appears to be highly superior to the procedure in (32).This applies even more in light of the fact that the computing times to solve a one-dimensional wave packet expansion problem are in any case very short.
In a two-or three-dimensional wave packet expansion problem the number of grid points is so large that a quadratic numerical effort is unacceptably high.However, we will see shortly that the numerical effort of the multi-dimensional discrete Green's function approximation is not quadratic in the number of grid points but in fact much lower.

Exploiting the separability of the Green's function
The Green's function formalism for the d-dimensional free Schrödinger equation reads [31] using the free-particle propagator for t > 0.
Analogously to the one-dimensional case we consider the initial and final computational domains respectively.Moreover, we define a discrete representation of the initial wave function on Ω 0 using the grid points x ,j = (j − J /2) x , x = L /J , j = 0, . . ., J − 1, = 1, . . ., d, where J 1 , . . ., J d ∈ N are assumed to be even integers.Likewise, we define a numerical approximation of ψ(•, t) on Ω.To this end, we introduce the grid points and the approximation where K 1 , . . ., K d ∈ N.With these definitions and in analogy to the one-dimensional problem the approximation of the exact solution ( 33) is given by The numerical effort of this approximation is in In a typical wave packet expansion problem the number of grid points J 1 , . . ., J d and K 1 , . . ., K d are of the same order of magnitude such that the numerical effort O(J 2 ) increases quadratically with the number of grid points Noting that the free-particle propagator (34) is the product of d one-dimensional free-particle propagators G (1) in ( 8), the integral (33) can also be written as and the approximation (35) becomes for = 1, . . ., d.In d steps, the initial wave function is transformed along the different spatial directions.By counting all elementary operations in (37a) we immediately see that the multidimensional discrete Green's function approximation can be evaluated in computational time.Consequently, under the assumptions in (36), the numerical effort is in ).In terms of the total number of grid points J = J 1 J 2 . . ., J d−1 J d we have J 1 ≈ J 1/d and hence the numerical effort is in O(dJ (d+1)/d ).Our main interest is in the spatial dimensions d = 1, 2, 3, in which case the numerical effort is in O(J 2 ), O(J 3/2 ) and O(J 4/3 ), respectively.
We refrain from deriving error estimates of the multi-dimensional approximation but would like to point out that the multi-dimensional Green's function method is based on the same integral approximation as the one-dimensional method.Therefore, we expect that the observations from the previous section are largely transferable to the multidimensional case.

22:
return ψ 23: end function The discrete Green's function approximation of the one-dimensional wave packet expansion problem is computed using a single dense matrix-vector multiplication for which highly optimized routines are available in practically any programming language.The multi-dimensional approximation (37) can be realized using a series of simple for-loops which, however, is very slow in many scripting languages like Matlab or Python.
An alternative approach is to combine all matrix-vector multiplications in the th step of the calculation in a single matrix-matrix multiplication.To do this, the tensor to be transformed must first be arranged in the form of an extremely large matrix using a reshape operation.This matrix is then multiplied from left by the matrix G in (37b).Finally, the result matrix is reshaped back into the form of a tensor.This process is described for the three-dimensional wave packet expansion problem in Algorithm 1. Effectively, the entire calculation is carried out using three large scale matrix-matrix multiplications, for which very efficient routines are available in any practically relevant programming language.The calculations in this work were implemented in Python using the packages Numpy [32] and PyTorch [33].Let us first consider the interference of three Gaussian wave packets for x 1 , x 2 , x 3 ∈ R and t ≥ 0 using  τ n = 2σ 2 n and n = 1, 2, 3.The wave function in (38) defines our initial condition at t = 0 and serves as a reference solution for t > 0. The parameters in the example are given by σ 1 = σ 2 = σ 3 = 0.4 and as the domains of the approximations to the expanded wave function ψ at the final times t = 2, t = 8 and t = 32, respectively.The results of the discrete Green's function approximation using J 1 = J 2 = J 3 = 256 and K 1 = K 2 = K 3 = 256 grid points are shown in Fig. 8 and Fig. 9 in form of false-color plots at x 2 ≡ 0 and x 3 ≡ 0. The first figure shows the density and the second figure shows the imaginary part of the wave function.Here, and in all false-color plots below, the densities and the imaginary parts of the three-dimensional wave functions have been normalized to their individual maximum (absolute) value.At the given resolution, the countless oscillations in the imaginary part of the wave function at t = 32 cause undesirable aliasing effects in the graphical representation.In order to make the high-frequency character of the expanded wave function visible, we have computed an additional approximation on the domain Ω = [30,50] × [30, 50] × [−10, 10] (at the same resolution ).This region is indicated by a yellow box in Fig. 8 (d).The corresponding imaginary part of the computed approximation can be seen in Fig. 9 (d).
The approximations depicted in Fig. 8 and Fig. 9 are practically indistinguishable from the exact solution (38).This is illustrated by Fig. 10 (a) in which we plot the relative error as a function of the parameter J 1 = 32, 34, . ., 256 using J 1 = J 2 = J 3 and Additionally, Fig. 10 (b) shows the runtime of the three-dimensional discrete Green's function method as a function of the parameter J 1 .We compare two implementations of Algorithm 1.In the first case, the algorithm is implemented on single core of a CPU8 (Numpy), while in the second case, the algorithm is implemented on a GPU9 (PyTorch).The GPU computations are incredibly fast.In fact, they are several orders of magnitude faster than the computations on the CPU.For example, using J 1 = 256 (J = J 1 J 2 J 3 = 16 777 216 grid points), the GPU calculation is 411 times faster than the corresponding calculation on the CPU.

Example: Interference of two ring-shaped wave packets
Finally, we consider the interference of two ring-shaped wave packets in three spatial dimensions.In this context, we first construct a ring-shaped solution to the free Schrödinger equation in three spatial dimensions.Here, the motion in the x 3 -direction is described by a simple Gaussian wave packet [24] ζ(x  using τ = 2σ 2 , while the motion in the (x 1 , x 2 )-plane is given by a less trivial ring-like function that still needs to be specified.According to (33) we have using polar coordinates r = x 2 1 + x 2 2 and θ = arctan(x 2 /x 1 ).From now on we we restrict ourselves to axially symmetric initial conditions, i.e., ξ(r , θ , 0) = ξ(r , 0) being independent of θ.Due to the rotational symmetry of the Schrödinger equation it remains depending on r only also at t > 0. By means of the well-known formula [34] 1 2π where J 0 (α) is the Bessel function of the first kind of the order 0, we obtain Next, we further restrict ourselves to initial conditions of the form where κ, q, and s are real parameters.If s > 0, the wave function at t = 0 has a ring-like structure with the maximum of its absolute value at r = s/2 κ −1 , where, without loss of generality, we assume κ > 0. The normalization 2π where Γ(z) is the gamma-function [34].We then apply the series expansion of the Bessel func-tion [34] J and perform, after substitution of (42) into (41), the integration over r term by term.Summation of the obtained series yields where is the confluent hypergeometric (Kummer's) function [34] and (a) k = Γ(a + k)/Γ(a).
Using the ring-like solution in (39) we define the wave function which defines our initial condition at t = 0 and serves as a reference solution for t > 0. In the example below we choose δ = 3, κ = 0.75, s = 10, q = 0.5 and σ = 0.85.Moreover, the initial wave function is evaluated on the domain respectively.Densities and imaginary parts of such approximations are shown in the form of false-color plots at x 2 ≡ 0 and x 3 ≡ 0 in Fig. 11 and Fig. 12, respectively.In the calculations, the initial condition was evaluated using J 1 = J 2 = 256 and J 3 = 128 grid points.The same number of grid points K 1 = K 2 = 256 and K 3 = 128 was used to approximate the expanded wave functions.Once again, it becomes evident how much the temporal development of the density ρ = |ψ| 2 differs from the temporal development of the imaginary part (ψ).While the density of the expanded wave function at t = 32 appears to be very smooth, the short wavelengths in the imaginary part of the initial condition are still visible in the imaginary part of the final approximation.
Even this example is solved with incredible accuracy, as shown in Fig. 13 (a) where the relative error is plotted as a function of the parameter J 3 = 16, 18, . . ., 128 with J 1 = J 2 = 2J 3 .The parameters K 1 = K 2 = 256 and K 3 = 128 remain constant throughout the whole simulation.For completeness, Fig. 13 (b) shows the computational runtimes which, due to the smaller number of grid points, are even shorter than in the previous example.

Conclusion
We have shown that the discrete Green's function approximation yields highly accurate numerical solutions to the free wave packet expansion problem at small numerical costs.This method will, for example, greatly simplify the simulation of expanding Bose-Einstein condensates after they have been released from the trap and at times when atomic interactions have become negligible.We are therefore convinced that this work represents an important progress in terms of modeling, planning and interpretation of experiments in matter-wave interferometry.
To prove Lemma 3 we first note that H n (−z) = (−1) n H n (z) for every z ∈ C and n ∈ N 0 .
Consequently, it remains to proof that f n is non-decreasing on [0, ∞).To this end, we will show that , where H n are the Hermite polynomials, R = q (L/2) + max(|a|, |b|) and q = 1/ √ 2t.Moreover, applying Leibniz's rule to f (x, x , t) = G (1) function ψ : R d × R → C in d spatial dimensions.The initial wave function ψ 0 is assumed to be a smooth function that is either compactly supported on

∥Figure 2 :
Figure 2: Real part (a) and imaginary part (b) of the asymptotic approximation ψ in (6) and the exact solution ψ in (3) at t = 8 for the initial wave function ψ0 in (2).Relative error (c) of the asymptotic approximation ψ(•, t) on the interval Ω = [−40, 40] as a function of time t.
(a) and (b) on the interval [−40, 40] for t = 8.Their relative difference on the same interval is shown in Fig. 2 (c) for t ∈ [0, 200].It is clearly visible that the relative error decreases not faster than O(t −1

Figure 4 :
Figure 4: Relative error of the discrete Green's function approximation applied to the initial wave packet ψ0 in (2) using two phase-shifted Gaussians of width σ.

Figure 5 :
Figure5: Convergence of the discrete Green's function approximation for the example of two phase shifted Gaussian wave packets using the parameters given in(22).Error bounds Γc(J) in (23) for c = c using c = 0.2 and = 1, . . ., 30.Also shown is the maximum difference between the discrete Green's function approximation ψ and the exact solution ψ evaluated at K = 1024 grid points in the final domain Ω = [a, b].

Figure 6 :
Figure 6: Wave packet expansion of a compactly supported initial wave function.The initial wave function ψ0 is given in (24) using the high-order polynomial u in (29).

Figure 8 :
Figure 8: Expansion and interference of three phase shifted Gaussian wave packets in three spatial dimensions at different times.The initial wave function is given in (38) using t = 0. (a) Density of the initial wave function.(b)-(d) Density of the approximations of the expanded wave function at t = 2, t = 8 and t = 32.

Figure 9 :
Figure 9: Expansion and interference of three phase shifted Gaussian wave packets in three spatial dimensions at different times.The initial wave function is given in (38) using t = 0. (a) Imaginary part of the initial wave function.(b)-(d) Imaginary part of the approximations of the expanded wave function at t = 2, t = 8 and t = 32.In (d) the approximation to the expanded wave function is computed on the region marked by a yellow square depicted in Fig. 8 (d).

Figure 10 :
Figure 10: Relative error (a) and run time (b) as a function of the parameter J1 = J2 = J3 for the example showing the expansion and interference of three phase-shifted Gaussian wave packets.

Figure 11 :
Figure 11: Expansion and interference of phase shifted ring-shaped wave packets in three spatial dimensions at different times.The initial wave function is given in (44) using t = 0. (a) Density of the initial wave function.(b)-(d) Density of the approximations of the expanded wave function at t = 2, t = 8 and t = 32.

2π 0 dθFigure 12 :
Figure 12: Expansion and interference of two phase shifted ring-shaped wave packets in three spatial dimensions at different times.The initial wave function is given in (44) using t = 0. (a) Imaginary part of the initial wave function.(b)-(d) Imaginary part of the approximations of the expanded wave function at t = 2, t = 8 and t = 32.

Figure 13 :
Figure 13: Relative error (a) and run time (b) as a function of the parameter J3 using J1 = J2 = 2J3 for the example showing the expansion and interference of two phase-shifted ring-shaped wave packets.