Spectra of Neutron Wave Functions in Earth's Gravitational Field

The time evolution of a quantum wave packet in the linear gravity potential is known as Quantum Bouncing Ball. The qBounce collaboration recently observed such a system by dropping wave packets of ultracold neutrons by a height of roughly 30 microns. In this article, space and momentum spectra as well as Wigner functions of the neutron wave functions in the gravitational field of the Earth are analyzed. We investigate the quantum states in the"preparation region", into which they transition after exiting a narrow double-mirror system and where we would expect to observe free fall and bounces in classical physics. For this, we start from the stationary solutions and eigenvalues of the Schr\"odinger equation in terms of Airy functions and their zeros. Subsequently, we examine space and momentum distributions as well as Wigner functions in phase space for pure and mixed quantum states. The eventual influence of Yukawa-like forces for small distances of several micrometers from the mirror is included through first order perturbation calculations. Those allow us to study the resulting modifications of space and momentum distributions, and phase space functions.


I. INTRODUCTION
A quantum wave packet bouncing on a hard surface under the influence of gravity has drawn some attention in the literature due to its departures from classical behaviour [1][2][3][4][5][6][7]. Other aspects of this quantum bouncer have also been studied to some extent. Of those, we would like to mention its chaotic behavior [8], the mathematical basis with orthonormal Airy eigenfunction solutions [9], the Wigner phase space as an interface of gravity and quantum mechanics [10], quantum revivals in a periodically driven gravitational cavity [11], and inertial and gravitational mass in quantum mechanics [12]. The development of sufficient ultracold neutron sources at the Institut Laue Langevin (ILL) in Grenoble and techniques to manipulate neutrons with high precision have made the simple quantum bouncer experimentally realizable. Demonstrations of quantum states in the gravitational potential of the Earth can be found in [13][14][15] and aspects from a more theoretical point of view in [16,17]. From the beginning these experiments were used to constrain hypothetical gravity-like interactions [18][19][20].
In this article, we will examine some details of the bounce of a neutron wave packet closely related to an experimental realization by the qBounce collaboration. More precisely, we will investigate the behavior of the momentum space wave packet solutions, the widths of the position and momentum space wave packets during the "bounce", and aspects of Yukawa-type interactions. Extensive use of the Wigner function formalism as a function of time is made as well.
The qBounce experiment has been performed at the UCN-beam position of the PF2 instrument at ILL, so far the 7th strongest source for ultracold neutrons with high continuous fluence, which is ideal for quantum bouncer realizations. It tests gravity at small distances with quantum interference techniques. The experimental tool is a gravitationally interacting quantum system -an ultracold neutron in the gravitational potential of the Earth -and a reflecting mirror above which the neutron is bound in well-defined quantum states. The collaboration is continuously developing a gravity resonance spectroscopy (GRS) [21][22][23][24] technique, which allows for a clear identification of the measured energy eigenstates states |1 −→ |2 , |1 −→ |3 , |1 −→ |4 , |2 −→ |3 , |2 −→ |4 , |2 −→ |5 , and most recently |1 −→ |6 . In this way, precisions are reached which enable us to search for hypothetical gravity-like interactions with relevance for cosmology. So far limits for axions [22], chameleon [23] and symmetron fields [24] have been placed.
For the purpose of this article, an important observable is the spatial density distribution of a free falling neutron above a reflecting mirror. A newly developed position-dependent neutron detector makes it possible to visualize the square of the Schrödinger wave function [25,26]. Detailed descriptions of these processes can be found in [27]. We now have a high-precision gravitational neutron spectrometer with available spatial resolution of 1.5 µm at our disposal.

II. SCHRÖDINGER EQUATION FOR QBOUNCE
The time-dependent Schrödinger equation for a neutron with mass m N in the gravitational field of the Earth with potential energy (g is the gravitational acceleration, z the distance above the mirror) whereĤ is the Hamiltonian containing V (z). The energy of the wave function ψ(z, t) is quantized in the potential V (z). Using the ansatz ψ n (z, t) = e −i En t ψ n (z) for a stationary state of energy E n with n = 1, 2, ..., we obtain the time-independent Schrödinger equation For negative values of z we have ψ(z) = 0 because the particles cannot enter the mirror surface. Therefore, the boundary condition for the solution of the differential equation is ψ(0) = 0. For this reason, it is supposed that the surface of the mirror has an infinite Fermi potential and the quantum wave does not enter the surface. This is of course just an approximation. At this point, it is appropriate to mention that the problem of two mirrors as well as the transition from an inertial frame (z 0 , t 0 ) to a non-inertial frame (z, t) has already been described in [29].
Searching for solutions ψ n (z), we multiply Eq. (4) with the factor and, using the substitutions [27] where E 0 is a characteristic gravitational energy scale, we obtain the differential equation Comparing this equation with the Airy equation we notice that the (non-normalized) eigenfunctions ψ n (ζ) can be expressed through the Airy function Ai(ζ) by moving the origin of coordinates to the n th zero point a n : ψ n (ζ) = Ai(ζ + a n )Θ(ζ) .
As described in [29] the normalized wave function ψ n (z, t) is given by The first zero point of the Airy function is located at a 1 ≈ −2.3381. This means, that E 1 = −a 1 E 0 ≈ 1.41 peV. Additional zero points are located along the negative axis, as can be seen in Fig. 1. They determine the energy eigenvalues E n according to Eq. (5). We list some of them below (z n = −z 0 a n ): The quantities z n are given, such that they can later be compared to Eq. (29). Concerning the calculation of the spectra and the Wigner function, the following formalism is developed using the example of the ground state. The wave function of the ground state ψ 1 (ζ) can be written as ψ 1 (ζ) = Ai(ζ + a 1 )Θ(ζ), where Θ(ζ) is Heaviside's step function, see Fig. 2a. This Heaviside step function is necessary in order to fulfill the boundary condition caused by the mirror whereupon the wave function has to be zero for negative ζ−values.
The spatial distribution is given by |ψ 1 (ζ)| 2 . This function can be taken from Fig. 1, orange curve, by imagining that the curve is shifted by a 1 to the positive ζ−axis. In order to attain the momentum space (variable k), the wave function ψ 1 (ζ) has to be Fourier transformed: where the two functions sin(ζk)Ai(ζ + a 1 ) dζ (12) have been defined. They are displayed in Fig. 3. The momentum spectrum is given by see Fig. 2b. We have to stress that k is actually a dimensionless variable and related to the physical momentum k p by B. Excited states Here we look at the excited states by discussing their momentum spectra for a few selected example values of n. The first excited state is characterized by the second zero point a 2 ≈ −4.08795 of the Airy function. Its wave function and momentum spectrum are depicted in Figs. 2a and 2b, respectively.
The third zero point of the Airy function is located at a 3 ≈ −5.52056 (2 nd excited state) and yields a momentum spectrum as given in Fig. 4. Besides, this figure also shows the results for the 3 rd and 9 th excited states with the corresponding fourth and tenth zero points a 4 ≈ −6.78671 and a 10 ≈ −12.8288 of the Airy function.
In Fig. 4 we can see that the number of oscillations before the onset of the asymptotic behavior of the momentum spectra increases with n. In case of n towards infinity, the amplitudes of the oscillations tend to zero and the momentum spectrum becomes a constant.

C. Presentation using Wigner function
The 2-dimensional Wigner function is an important tool in quantum optics. It allows for a simultaneous view into space and momentum regions. The Wigner function is real but can be positive and negative as well. In this respect, it is not a classical 2-dim distribution function. Therefore, it is often called a quasi-distribution function. Most remarkably is the property that an integration of a Wigner function over momentum gives the spatial probability, while integration over the spatial coordinate gives the momentum probability. These two marginal distribution functions (spatial and momentum distributions) are, at least in principle, experimentally accessible. The Wigner function formalism has already been applied within the framework of investigations of the gravitational potential of the Earth [12]. In addition, we would like to point to an article, in which the interface of gravity and quantum mechanics has been discussed with the Wigner phase space distribution function [30].
The Wigner function in Eq. (16) can be rewritten as: Using Eq. (18), we can find the momentum distribution by integrating over ζ: Evaluating this expression numerically for the states considered in Figs. 2b, we obtain the same results as in these figures. Consequently, the Wigner function provides us with a second option to calculate |F (k, a n )| 2 . The 2-dim Wigner function Eq. (18) of the ground state W (ζ, k, a 1 ) is plotted in Fig. 5a. It is almost everywhere positive. There are only very small and hardly visible negative regions. The spatial probability density has been experimentally verified for ultracold neutrons in [23]. Here we suggest that the momentum probability distributions could be measured in a similar fashion. For example, if the ground state population amounts to 70% (p 1 = 0.7), the first excited state amounts to 30% (p 2 = 0.3), and no other excited states are populated, then extracting the total probability distributions from Fig. 2b is straightforward. This is because the relative contributions can be extracted from the figures: n p n |F (k, a n )| 2 , see Fig. 6. The orange line in Fig. 6 presents an example for p 1 = p 2 = 0.5. In this case, the first excited state, represented by |F (k, a 2 )| 2 and Fig. 2b, is clearly visible. It should be mentioned that this procedure corresponds to an incoherent superposition. Therefore, the proposed procedure is not exact. We have to take the time dependence of the wave function, see Eq. (3), into account. Consequently, we will now go back to the time-dependent space distribution |ψ(ζ, t)| 2 using the following ansatz of coherent superposition where ψ 1 (ζ) = Ai(ζ +a 1 ) Θ(ζ) and ψ 2 (ζ) = Ai(ζ +a 2 ) Θ(ζ) are real functions, and we ignore a potential phase between both terms for simplicity. The time-dependent position probability of this superposition state is therefore This function oscillates with time t. If p 1 = 1 and p 2 = 0, we recover |ψ 1 (ζ)| 2 , which is the square of the function depicted in blue in Fig. 2a. Furthermore, we have This means, that |ψ s (ζ, t)| 2 oscillates in the time-range of milliseconds, as can be seen in the example shown in Fig. 7a. The first small peak at t = 0 appears at ζ ≈ 2.5 and not at ζ ≈ 1.3, i.e., the maximum of ψ 1 (ζ), see Fig. 2a, because the interference term in Eq. (21) contains Ai(ζ + a 2 ), which is negative in the region between ζ = 0 and ζ ≈ 1.7, see Fig. 2a. However, the large maximum at t ≈ 0.002 s appears for ζ ≈ 1.3 due to ψ 1 (ζ) and p 1 = 0.7. We assume that such properties of time-dependent position probabilities for superpositions of ground and excited states can be measured.
We have to point out that in case of a non-time-resolving measurement we have to integrate time t in Eq. (21) over one time period, e.g., from 0 to 2π /(E 1 − E 2 ). In this case, the interference term disappears and we obtain only the first two terms Fig. 8 depicts an example. This spatial probability density is very similar to the results of an experiment using a track detector [23].
The Fourier transform of ψ s (ζ) in Eq. (20) reads Next, we want to calculate |F s (k, t)| 2 . We obtain  Using Eq. (11) in order to decompose β * , we find Since The time dependence of this momentum spectrum could also be measured experimentally. Fig. 7b gives an example of Eq. (26) with p 1 = 0.7 and p 2 = 0.3. If we chose p 1 = 1 and p 2 = 0, or p 1 = 0 and p 2 = 1, we would instead recover Fig. 2b. In the case of a non-time-resolving measurement, we have to integrate time t in Eq. (26) over one time period. Hence, the interference term disappears and we obtain only the first two terms p 1 |F (k, a 1 )| 2 + p 2 |F (k, a 2 )| 2 , which, using p 1 = 0.7 and p 2 = 0.3, is the function in blue in Fig. 6.
Using Eqs. (16) and (20), the Wigner function of the coherent superposition can be found to be The single Wigner functions W (ζ, k, a 1 ) and W (ζ, k, a 2 ) have already been depicted in Figs. 5a and 5b, respectively, while the integral on the right-hand side of Eq. (27), the interference term, can be evaluated similarly to Eq. (19). If the time t cannot be resolved experimentally, we have to integrate over one time period, which causes the interference term to disappear. The resulting Wigner function for an example population is given in Fig. 9. L from article [32] equation (10): Here the prime denotes derivatives with respect to the argument, i.e., z 0 d/dz, and Ai(x) and Bi(x) are the two independent solutions of Airy's equation. Due to the experimental setting, we assume that the wave function in region I has support only on [0, L] × (−∞, 0], but we keep this assumption implicit for notational convenience. Below we present a selection of possible numerical values for the parameters used in Eq. (28): Here we haveĒ m =z m m N g. The energy spectrumĒ m is obtained by the conditions that the wave functions vanish at the lower as well as upper mirror surface, i.e., ψ A. Fourier transformation of the wave function The Fourier transformation of the wave function in Eq. (29) is given by where C m (t) = 1 √ 2π √ z0Nm e − i Ē mt and we made use of the restrictions on the support of the wave function that we mentioned earlier. Since z has the dimension of a length, here the variable k must have the dimension of an inverse length and is related to the physical momentum k p by We define the following stationary quantities: such that and the spectral function is given through with |C m | 2 = 1 2πz0N 2 m . Notice that this momentum distribution is stationary. It is depicted for the cases m = 1, m = 2, and m = 3 in Figs. 11a, 11b, and 11c, respectively. These figures also explicitly show the respective |C m | 2 α 2 s (k, m) and |C m | 2 α 2 c (k, m). The wave function for two mirrors ψ In a similar way, the spectral function |F We can easily show that this normalization is indeed valid here by using the first equality in Eq. (30): and A(z) and B(z) are the limits of integration that appear due to the restrictions on the support of ψ (0) m . It can directly be seen that D(z, z ) = D(z, −z ).
We will now determine the limits of integration. For this, we consider that z is limited between 0 and L, leading us to the following 4 equations which represent 4 straight lines in the (z, z )-diagram and generate a rhombus. The values inside of this rhombus are the allowed values for integration. Therfore, we conclude Since A(z) = −B(z), the Wigner function can finally be written as It is depicted in Figs. 12a, 12b, and 12c for m = 1, m = 2, and m = 3, respectively. As can be seen from Fig. 12a, the Wigner function W  In summary, in this chapter, we considered the eigenstates of the neutron wave function in a double mirror system, calculated the spectral functions and evaluated the corresponding Wigner functions. The latter enabled us to look into the complete phase space in order to study momentum and position simultaneously.

IV. "FREE FALL" OF WAVE FUNCTION AFTER DOUBLE MIRROR (REGION II)
In this chapter we consider the "free fall" of a wave function which exits a double mirror system (we denote this region by I, see Fig. 10). The wave function reaches a second region II, where it falls down a height h = 27 µm on a subsequent static mirror located below the double mirror system. This case has been investigated theoretically in [29]. The wave function in region I has been given in Eq. (28). Since we are now also going to consider region II, and for convenience, we apply the coordinate shift z → z − h to the result from Eq. (28), such that we have where we introduced the notationC m := 1 √ z0Nm , and the wave function has support only on [h, L + h] × (−∞, 0]. In region II the wave function takes on the following form and includes coefficients (47) Fig. 13 shows some of the first coefficients D n,m for the cases m = 1 and m = 2. From there it can be seen that for n > 12 the coefficients D n,m become very small and can therefore be neglected. Now we can check whether ψ 1,I (z, t = 0) and ψ 1,II (z, t = 0) fulfill the condition in Eq. (46). We do this for the ground state m = 1. For this, we consider only a finite number of coefficients D n,1 . The results including n up to 12 and 15 are presented in Figs. 14a and 14b, respectively. There the abscissa represents the z−coordinate in µm shifted to the point of origin of the coordinate system. Therefore, in region I there are values of z from 0 to L. The agreement between ψ 1,I (z, t = 0) and the plotted approximation of ψ 1,II (z, t = 0) is very good except near certain regions, e.g., around z = 0. These small differences supposedly are due to using only a finite number of D n,1 -coefficients. We expect the differences to become smaller when more D n,1 -coefficients are considered.

A. Spatial distribution (SD) in "free-fall" region
The spatial distribution in region II follows from Eq. (44):   We define such that the spatial distribution in region II reads Whenever performing numerical calculations we will only consider the sums in Eq. (49) from n = 1 to 15 due to the smallness of later coefficients D n,m . In Figs. 15a and 15b the spatial distributions (SD) |ψ 1,II (z, t)| 2 and |ψ 2,II (z, t)| 2 , respectively, are plotted as a function of the z−coordinate (in µm) and time t (in s). For t = 0 the ground and the first excited state are visible. This means, that between z = 0 and z = h = 27 µm the SD are zero. Between z = h and z = h + L = 55 µm the SD have the shape of the ground state, see Fig. 15a, or the first excited state, Figs. 15b. For z > 55 µm the SD vanish again. While t evolves, the wave function is reflected from the mirror in region II multiple times. Since the frequencies E n / vary during this process, which leads to varying superpositions of waves, more complicated SD pictures result at times after t = 0 . Figs. 16a and 16b show the cosine and sine distribution functions |G c 1 (z, t)| 2 and |G s 1 (z, t)| 2 of the ground state. The sum of these functions yields Fig. 15a. It is interesting to consider these parts of |ψ 1,II (z, t)| 2 separately since they could be important for interpreting experimental results.
(b) Spatial distribution |ψ2,II(z, t)| 2 ; The first excited state (m = 2) appears at t = 0. So far, we have considered only a particular state (ground state or first excited state) which enters region II. It would be interesting to consider, for example, a mixture of these two types of quantum states arriving at the step before entering region II. This case will now be investigated. In the following, we distinguish between coherent mixtures and incoherent mixtures. We consider mixtures of the ground state with the first excited state.

Coherent mixtures
In the beginning, we look at the following coherent superposition state: where p 1 and p 2 are probabilities, and we ignore a potential phase between both terms for simplicity. Note that ∞ 0 |ψ 1+2,II (z, t)| 2 dz = 1 because of the orthonormality of the eigenfunctions: ∞ 0 ψ 1 ψ 2 dz = 0. The next step is to calculate the spatial distribution of the coherent superposition |ψ 1+2,II (z, t)| 2 . This can be

Incoherent mixtures
Incoherent mixtures can be described by the following formula: In Figs. 18a and 18b two examples are presented, for which incoherent mixtures can be observed, namely for p 1 = 0.7 and p 2 = 0.3, and p 1 = 0.5 and p 2 = 0.5, respectively. In Fig. 18a, for example, at t = 0 the shape of the SD is qualitatively the same as in Fig. 8. Altogether, the oscillations are less distinct than in the previous case. Figs. 19a and 19b show two examples, in which |ψ incoh mix,II (z, t)| 2 is separated into a cosine and a sine part. The formulas are included in the corresponding figure captions, and the sum of both of these plots gives the result in Fig. 18a.  The momentum distribution in region II is calculated through a Fourier transformation: This can be expressed as In Figs. 20a and 20b we present two examples of |F m,II (k, t)| 2 for m = 1 and m = 2. Fig. 20a shows the time dependence of the ground state momentum distribution in region II. At t = 0 the function |F 1,II (k, 0)| 2 is the Fourier transform of the approximated wave function of Fig. 14b squared and has a maximum exactly at k = 0. During the evolution of time t the maxima are displaced to larger values of k. This causes an inclined periodic (t − k)-pattern of the momentum distribution. For the first excited state (m = 2) the pattern is shifted once more, as can be seen in Fig. 20b. At t = 0 we can observe a double peak in accordance with the first excited state at the boundary between regions I and II. As expected, there is a distinct minimum at k = 0. The Wigner function W m,II (z, k, t) in region II is given in terms of the corresponding wave function ψ m,II (z, t) from Eq. (44): Proceeding similarly to how we did in Eq. (40), we find the limits of integration to fulfill −2z ≤ z ≤ 2z. Moreover, we can separate the wave function into two parts using Eq. (49): It is interesting to include a small perturbation to the potential in the Schrödinger equation in order to simulate a small variation of the gravitational field near the mirror. The basic idea behind this proposal is to modify gravity at small distances and determine the limits on non-Newtonian gravitation below 10 µm.

A. First order perturbation calculation and a new wave function
When a neutron with mass m N approaches the mirror, the mass of this extended source might modify the gravitational acceleration of Earth g due to a potentially present non-Newtonian force with range δ, see Ref. [21]. This modification would lead to an additional, Yukawa-type interaction The parameter W 0 is a positive or negative constant with the dimension of an energy and δ is called Yukawa-distance over which the corresponding force acts. δ is measured in units of z 0 . The stationary Schrödinger equation (4) reads now Ψ Yu n (z) and n are the corresponding wave functions and energy eigenvalues. The quantity W (z) has to be a small correction to m N gz. In the basis of normalized Airy-functions given in Eq. (9) we obtain at first order: where with J n ,n : Note that ψ n |ψ n = 1, ψ n |ψ (1) n = 0, ψ n |Ψ Yu n = 1 and Ψ Yu n |Ψ Yu n ≈ 1. In Fig. 22a the total potential V (z) = m N gz + W (z) is drawn for W 0 = 0 and for an additional attractive potential with the strength W 0 = −1 peV. Since the total force is therefore given by K(z) = − ∇V (z) = [−mg + W0 δ e −z/δ ] e z , we can see that the Yukawa-potential leads to an additional force. Fig. 22b shows n as a function of −W 0 . The eigenvalues n decrease with increasing W 0 .
such that An interesting way of comparing space distributions with and without Yukawa interaction is the following quantity: Here |ψ m,II (z, t)| 2 is given in Eq. (50) or, equivalently, by Eq. (76) when setting W 0 = 0. In Fig. 23a we assume W 0 = −1 peV and draw the corresponding ∆ Yu (z, t). For comparison, we also depicted |ψ Yu m,II (z, t)| 2 in Fig. 23b. For t = 0 this function is exactly the same as in Fig. 15a. However, for t > 0, small differences to the function shown in Fig. 15a can be noticed.
Setting W 0 = 0 in Eq. (79) recovers Eq. (57), which gives the expression for |F m,II (k, t)| 2 . An interesting way of comparing the momentum distributions with and without Yukawa interaction is the following quantity: In Fig. 24a this difference of the momentum distributions with and without Yukawa interaction is depicted. At the jump discontinuities t ≈ 3 ms and t ≈ 9 ms distinct differences are visible. At these times the wave function is reflected at the mirror where the Yukawa potential is the strongest. This can also be observed in Fig. 24b, in which |F Yu m,II (k, t)| 2 is shown. At these points in time |F Yu m,II (k, t)| 2 exhibits distinct maxima compared to the momentum distribution without Yukawa interaction depicted in Fig. 20a. In this theoretical treatise the wave function of the qBounce experiment has been investigated in detail. The gravitational field of the Earth constitutes the potential used in the Schrödinger equation. This yields solutions for the wave function which correspond to the Airy function. Since the wave function has to vanish at the mirror surface, a ground state and excited states evolve. These states have been analyzed with respect to spatial and momentum distributions. For this purpose, the Wigner function has also been used. It was shown that the distribution spectra of the ground and excited states exhibit the anticipated properties, which are reproduced in the marginal distribution functions. Furthermore, the time dependence of a mixture of the ground state and the first excited state has been considered. The qBounce-problem in which the neutron wave is enclosed between 2 mirrors has been analyzed as well. Finally, the case where the wave function exits the double mirror system and freely falls on a subsequent mirror has been considered. The purpose of these calculations was to motivate measurements both in real space and in momentum space for comparison between experimental findings and theoretical results. Finally and in addition, we made an attempt at a first order perturbation calculation in order to describe a very small change in the potential near the mirror due to a Yukawa-like coupling. Already from this very simplified calculation we predicted differences in the spatial and momentum distributions between cases with and without a Yukawa-like interaction. However, in order to make a statement about a realistic experimental situation, the probability distribution of neutrons at the transition from region I to region II should be known in detail when also taking into account the Yukawa-like interaction in region I. Though, this is beyond the scope of this article since a much more intricate computation would be required.