Scale-free relaxation of a wave packet in a quantum well with power-law tails

We propose a setup for which a power-law decay is predicted to be observable for generic and realistic conditions. The system we study is very simple: A quantum wave packet initially prepared in a potential well with (i) tails asymptotically decaying like ~ x^{-2} and (ii) an eigenvalues spectrum that shows a continuous part attached to the ground or equilibrium state. We analytically derive the asymptotic decay law from the spectral properties for generic, confined initial states. Our findings are supported by realistic numerical simulations for state-of-the-art expansion experiments with cold atoms.


Introduction
The temporal evolution of an initially localized, quantum mechanical wave packet is of fundamental importance for understanding the difference between classical and quantum transport, and goes back to Schrödinger's early attempts to construct wave groups which behave like classical particles [1]. With increasing complexity of the potential landscape, which can be due to its topography as well its topology, a panoply of surprising transport phenomena emerges, from ballistic over sub-, super-or just diffusive [2,3] to various types of localized transport [4][5][6][7]. These phenomena have their respective effective descriptions, e.g. in semiclassical [8], mesoscopic [9] or statistical [10] terms, and manifest in distinct experimental settings, from light-matter interaction [11] over fundamental (quantum) optics [12,13] to ultracold matter [14,15], quantum walks [16] and biochemistry [17]. On the fundamental level, however, all transport properties are hardwired in the spectral properties of the underlying quantum system, and control of the latter implies control of the former. Given the stunning experimental control over potential landscapes in variable dimensions, as achieved e.g. in cold matter science over the last decade [15,[18][19][20][21], we can contemplate exploring the above diversity of quantum transport phenomena for optimal control, by tuning the decisive spectral properties.
In this paper, we study the particularly simple, although paradigmatic case of a tunneling escape from a one-dimensional (1D) potential well. The proper choice of the confining potential allows one to induce algebraic rather than exponential decay, for generic, confined initial states. As described in [22,23], the crucial ingredient will be an asymptotically scaleinvariant potential decaying to zero as ∼x −2 , which occurs naturally in systems with dipolar interactions [24] or anomalous molecular binding potentials [25], and is often associated with quite counterintuitive effects (see, e.g., [22,[26][27][28][29]). We thus provide in section 2 a particularly transparent example of a scale-free relaxation process, where the equilibrium state is reached only at asymptotically long times [3,30]. The scenario considered here is shown to be robust against unavoidable experimental modifications of the idealized theoretical scenario we depart from. In section 3 we argue that algebraic, rather than exponential decay, should be observable in modern experiments with expanding Bose-Einstein condensates [15,20] in engineered potentials [18,19,21]. Section 4 concludes the paper.

Power-law relaxation
We start our analysis by considering the two-parameter family of 1D, single-particle potentials [27] We use dimensionless units in the sequel, measuring actions in units ofh, and giving the particle unit mass. The exact analytical solution of the Schrödinger equation for the potential V S is available [27,31]. We will, however, see that, given its asymptotic scaling, the precise form of the confining well is not crucial for the predicted power-law decay of the survival probability. The potential of equation (1) has the same asymptotic behavior as in [22,23]; however, differently from the one considered in [23] it extends over the entire real axis, implying the existence of an equilibrium or ground state.

Theoretical predictions
Equation (1) ensures that the ground state eigenfunction ψ 0 associated with the eigenvalue E = 0 is continuous with a continuous first derivative. For |x| > L one finds that [22,27], which ensures the square integrability for α > 1. For |x| L, the eigenfunction of the ground state is ψ 0 = B 0 cos( √ V 0 x). The constants A 0 and B 0 are set by imposing that ψ 0 is normalized to unity and continuous in x = ±L. The continuity of the first derivative in x = ±L implies the following relation between the three original parameters of the potential: Besides the discrete eigenvalue E = 0, the spectrum of V S has a continuous component in the range E > 0, which is attached to the null eigenvalue. Since the potential is even, we can have odd and even solutions to the Schrödinger equation. For |x| > L the eigenfunction ψ (odd) E is a linear combination of the Bessel functions ψ (odd) The coefficients a E , b E and d E are fixed by imposing that ψ (odd) E and its first derivative are continuous in x = L and that ψ (odd) E are orthonormalized with a δ-function of the energy: dx Similar conditions apply to the even solutions. Further details can be found in the appendix.
Let us now consider the particle prepared in an initial state 0 (x) = (x, 0), given by a linear combination of eigenfunctions {ψ 0 (x), ψ E (x)}, with real coefficients {a 0 , a E }. Its time evolution will be given by where we explicitly assume that a 0 = 0. For our purpose the distribution of participating energies in the initial state should be continuously connected to zero and be sufficiently broad, as we will discuss below in more detail. In contrast to [23], we are explicitly interested in the survival probability P(t) that the particle remains confined within the well, i.e.
By substituting equation (3), we obtain where For the 1D potential of equation (1), the integrals I 1 , I 2 and I 3 can be solved analytically, while this is, in general, not the case for the time-dependent integrals C 2 (t) and C 3 (t). However, we are interested only in the long-time behavior of P(t), and can therefore take advantage of tools from asymptotic theory of Laplace and Fourier transforms [32,33], to obtain the asymptotic form of P(t). The functions C 2 (t) and C 3 (t) can be written as Fourier integrals of the type with j = 2, 3. By considering the explicit form of the eigenfunctions ψ 0 , ψ E (see below), one gets for small values of E: with β = (α − 3)/2. Hereafter, we will only consider initial conditions 0 (x) for which a 0 = 0 and the spectral decomposition involves continuum energies that extend down to E = 0.
The explicit values of C 2 and C 3 clearly depend on the initial condition. Hereafter, we will give some results for a wide range of initial conditions. For initial states that take non-vanishing values only in the well [−L , L], a E ≈ E β in the limit E → 0 [27]. Equation (8), furthermore, allows us to infer for small values of E, where K (α) and H (α) are prefactors that depend also on the initial condition. If the initial state 0 (x) has non-vanishing values outside the well region [−L , L], the above results still hold, provided that, for large |x| L, 0 (x) exhibits power law ∼|x| −a , Gaussian ∼ exp(−ax 2 ), or stretched exponential ∼ exp(−ax b ) decay. The leading contribution to the asymptotic long-time behavior of the survival probability is then found to be associated with the C 2 (t) term above, and reads Through the dependence of K (α) on the spectral expansion of 0 (x), via f 2 or C 2 , the asymptotic decay of P(t) also depends on the initial condition. Thus, as anticipated in the  (7), for the potential (1) and different initial conditions: introduction, we have shown that the potential of equation (1) induces an algebraic decay in the survival probability P(t), for generic, confined initial states. In contrast, for an initial state 0 (x) with a E = 0 for E < E c (for some E c > 0), one can show that the survival probability has an exponential cut-off e −E c t at large t, which is a manifestation of the power-law decay arising from continuous spectral components arbitrarily close to the ground state energy. Figure 1 shows P 2 for the initial condition 0 (x) = δ(x). For some values of α, e.g. in α = 1 and for α = 3, P 2 vanishes. The oscillatory behavior of P 2 stems from the trigonometric functions in (11), while the fact that P 2 becomes negligible for α 7 is induced by K (α), and is therefore related to the initial conditions as used in figure 1. The oscillatory behavior of P 2 may result in a non-monotonic decay of the survival probability (an example will be shown in figure 3). In fact, a negative value of P 2 implies that the asymptotic value |a 0 | 2 I 1 is approached from below and therefore the probability derivative must change in sign. This is not peculiar to the system we are studying here. In fact, a similar behavior is also observed starting from an usual square well quantum potential V S (x) = −V 0 when |x| L and V S (x) = 0 when |x| > L, with 0 (x) = δ(x). For vanishing P 2 , one has to consider the next leading time-dependent term in the asymptotic expansion of C 2 and C 3 .
In [23], according to equation (35) therein, a specific initial state was chosen whose expansion coefficients a E show a power-law behavior, although with an exponent different from ours, for E → 0. This implies a power law for the survival probability which is different from our result of equation (10). The main difference between our potential of equation (1) and the 1D potential considered in [23] is that the latter is defined only on the positive real axis, while our potential extends over the whole real axis. Consequently, in [23], there does not exist a bound state at E = 0, leading to a decay without a lower bound, or in other words without control.
In contrast, our equation (10) explicitly takes into account a 0 = 0. Another difference as compared to [23] is that we explicitly study the survival probability inside the quantum well, as defined in equation (4) and in view of our proposed experiment in section 3.3, and not the fidelity function considered in [23].

The role of the potential tails
In this section, we will show that the mere absence of a gap between the null eigenvalue and the continuum part of the spectrum is a necessary but not sufficient condition to observe a powerlaw decaying survival probability.
In general, any Schrödinger potential that asymptotically decays to zero like 1/x µ has a continuum part of the spectrum attached to the null eigenvalue. Let us now, therefore, consider the more general case of a Schrödinger potential where L, V 0 , V 1 and µ are real positive constants. Let us call {ψ (µ) 0 , ψ (µ) E } the eigenfunctions associated with such potential. In this case, we can only obtain the eigenfunction associated with the null eigenvalue: The eigenfunctions relative to the continuum part of the spectrum are not known. The parameters L, V 0 , V 1 and µ can be chosen in such a way that the spectrum contains a single discrete eigenvalue E 0 = 0 and a continuous part for E > 0. As a result, the parameters L, V 0 , V 1 and µ are not independent. In fact, the continuity of ∂ x ψ 0 in x = L provides a relation between them: In the following, we consider L, V 0 and µ as independent parameters and will obtain V 1 by numerically solving the above equation.
We will prove below that the survival probability of the above process of equation (12) is not, in general, power-law decaying. In fact, we will show that the C 2 term may eventually decay like a power law for large time values only for very specific initial conditions. Let us make two preliminary observations. Firstly, it is worth mentioning that by using the orthogonality relation one can prove the following identity: for I 2 as defined in equation (6). Secondly, we notice that in order to have C 2 (τ ) ≈ 1/τ a+1 one must have that a E I 2 (E) ≈ E a for small energy values.
Let us now consider the large x behavior of For large values of x the potential vanishes; thus ψ (µ) x . Therefore, by using equations (16) and (12.01) in chapter 3 of [32], the ansatz This result implies that ψ (µ) 0 (x) should decay according to a power law, given equation (16). However, by using equation (13), one has that showing that I 2 (E) cannot behave like a power law for small values of E when µ < 2. It is worth mentioning that when µ = 2, the above ansatz holds true with a = (α − 5)/4. In this case, one would obtain I (x) ≈ x 1−α/2 , which is in agreement with the fact that ψ 0 (x) = x −α/2 for large x values. We have, therefore, shown that I 2 (E) is not growing like a power law for small energy values. Therefore, in order to have a power-law decay in C 2 (t), one has to engineer appropriate initial conditions such that a E I 2 (E) behaves like a power law for small energy values. In conclusion, the absence of an upper bound for the time scale is a necessary but not sufficient condition in order to observe a power-law decaying auto-correlation function. This also implies that the decay of the survival probability is faster for µ < 2 than for the case with µ = 2.

Numerical simulations of the model
Because of its discontinuities, it is difficult to reproduce V S (x) experimentally or in numerical simulations. We therefore consider the smoothed version V sm (x) defined as [26] V sm (x) = α 4 in our subsequent numerical tests. Also this potential has a continuous spectrum attached to the null eigenvalue, and decays like x −2 at |x| → ∞. We are interested in the relaxation properties of wave packets prepared inside the potential well (see figure 2). We define the survival probability over the region bounded by the potential maxima at x c : Our observable is the approach of the survival probability to its asymptotic constant value P ∞ = |a 0 | 2 I 1 described by equation (10), i.e. the quantity Even if we simulate the evolution of the initial wave packet only in one dimension, the numerical computations take quite a long time for several reasons. First of all, we need to propagate sufficiently far out into the tails of the potential, while-at the same time-the expansion of the decaying parts is rather fast due to the high-energy components of the initial state. This is also the reason why simple absorption methods at the numerical boundary do not work very well since those are typically adapted to absorb just a small window of energies with sufficient precision. Moreover, in order to estimate the asymptotic probability P ∞ , which is not analytically known for the smooth potential V sm , we must propagate considerably longer than shown in the following figures, for which we approximate P ∞ by P(t max ), with t max = 140-200. For the actual observation of the predicted asymptotic behavior, the precise form of the initial state prepared in the well is not crucial, provided that it is given by a coherent superposition of energy eigenstates from a continuous energy range including the ground state energy. We consider a smooth initial state where N is a normalization constant. The problem defined by the potential of equation (19) is not exactly solvable unfortunately. Therefore, the expansion coefficients {a 0 , a E } are not known. This prevents us from giving an analytical solution for the survival probability, given the above initial state. We are, therefore, forced to perform numerical simulations, as detailed below. However, the fact that the expansion coefficients are not known makes it also difficult to estimate the initial population of the states in the continuum part of the spectrum. We can, however, estimate the mean energyĒ of the initial state of equation (22). By using the parameters of figure 2, we get, for instance,Ē ≈ 2.6. The initial state (22), shown in figure 2, is numerically propagated in real time using an implicit norm-preserving Crank-Nicolson integration scheme [34], which also controls the boundary conditions very well [35]. Figure 3 reports the results of the numerical simulations. One can observe a clear algebraic decay as predicted by equation (7), for various values of α which induce different algebraic decay exponents. The earlier the power-law decay emerges, the larger the weight factor P 2 in (7). The asymptotic value P ∞ is approached from above (α = 2.8) or from below (α = 3.2, 4), depending on this factor's sign (compare figure 1). In the latter case, a change of sign of P(t) − P ∞ at finite times induces the discontinuities in the  (7).
first derivative on the logarithmic scale of figure 3. For a fixed α, the power-law scaling regime may be further enhanced by adapting the precise form of the initial state (e.g. by the above parameter σ ) and also by changing the second parameter γ of the potential of equation (19). Both enter our approximate formula for the prefactor P 2 , see equation (11) and figure 1. The figure explicitly shows how the decay of the survival probability may not necessarily be monotonic as, for example, in the green and blue curves. The effect is magnified by the logarithmic scale. As mentioned in section 2.1, this is connected with the prefactor P 2 , which inverts its sign, as shown in figure 1. We do not have a precise physical explanation for this oscillatory behavior. Indeed, the time evolution of ψ E is highly non-trivial. This induces an alternation between (i) phases when the bulk of the state is outside the well, and therefore we observe a depletion of the survival probability, and (ii) phases when there is a re-entrance of the state within the well, which causes a partial restoration of the survival probability. In fact, no matter what its mean energy is, 0 is the linear superposition of eigenstates, some of which have very low energies, in particular smaller than V S (x c ). Exactly these low-energy components give rise to the observed power law in the survival probability after the transient, as we have shown in section 2.1. The oscillatory transient behavior is also present for a simple squarewell potential, which tells us that the transient should mainly be caused by the high-energy components contributing to the dynamics while the initial state relaxes in the well.

The role of the tail in the potential
In the following, we will show that the predicted power-law decay is robust with respect to realistic experimental situations. The only parameter which needs to be controlled very well is the exponent of the potential tails, i.e. µ in the asymptotics of the potential scaling as ∼x −µ .   (24): (P(120) − P(t))/ t versus 100 − t = t − 20, with t ≡ 120 − t. Hence, we observe a relaxation with a fast decreasing slope as t → 0, in the special case of µ = 2. In all the other cases, the slope remains almost constant over the shown times, implying a much faster decay with time.
We use the form which for µ = 2 recovers V sm (x) from equation (19). For µ = 2, the dynamics of the system will be qualitatively very different, as we show in figure 4 directly for the survival probability. Instead of an asymptotically slow saturation toward the ground state within the potential, we observe a fast decay-with approximately constant slopes-for all values µ < 2 along the same time scales as in figure 3.
As described above around equation (21), it is computationally hard to estimate the saturation value of the survival probabilities, since the ground states are not known analytically and the data for µ = 2 imply that a reliable numerical estimate will be possible only after much longer propagations than shown in the figures. This may cause the wrong impression suggested by the main panels of figure 4: the ground state is approached in shorter absolute time for µ = 2 as compared to the other values of µ, yet what counts is the rate of approach, which is always smaller for µ = 2. This is what we intend as 'slower' decay in the latter case and is the reason for plotting the slopes defined by (P(120) − P(t)) / t (P(120) − P (20)) /100 (24) in the insets, with t ≡ 120 − t. In order to obtain comparable values, we divided by the extremal values, see the denominator of (24), which all occur at the minimal time t = 20 (chosen after non-universal transients at still smaller times). Only in the special case of µ = 2, the relaxation to the steady-state value becomes much slower with time, see the fast decreasing slopes in the insets (black solid line). The shown numerical results confirm the theoretical discussion above in section 2.2, predicting in essence a faster approach to the ground or equilibrium state for µ < 2. Consequently, a possible experiment should control the exponent to be equal to 2 at least at 2-3 significant digits, in order to observe the power-law decay discussed in sections 2.1 and 3.1.

Experimental realization of the model
Let us conclude with an experimental protocol to test our prediction (10). We have in mind a Bose-Einstein condensate prepared in an optical trap, which is then exposed to a potential of the form of equation (19) while the trap is relaxed. The potential may hereby be created optically, for instance, with a fast-moving laser beam [18] or by holographic techniques [36]. A Gaussian initial state of the form is initially prepared in a harmonic trap with the characteristic oscillator length σ trap, i . N is a normalization constant. Then the smooth potential (19) is switched on, while we switch off or relax the trap abruptly to a shallow confinement characterized by the harmonic oscillator length σ trap, f σ trap, i . In the case of figure 5, we consider values from σ trap, f = 100 to 800, and σ trap, i = 0.5. In the absence of a trap (blue thick lines), the Gaussian initial state exhibits a behavior similar to that of the smooth water-bag state from above, equation (22). A shallow trap with large σ trap, f manifests at long times, by an exponential cut-off of P S (t), as shown in figure 5, since the potential's asymptotics are changed by the trap. This induces a spectral gap between the ground state energy and the continuum component, while it is specifically the absence of the gap which is responsible for the algebraic decay, as discussed in section 2.1. The steeper the confining trap potential, the larger the trap-induced spectral gap, and the shorter the time interval over which an algebraic decay can be observed, before the exponential cut-off. Indeed, we observe such a continuous degradation of the asymptotic law with exponent (α − 1)/2 = 1.5 when making the additional confinement steeper. However, the power-law trend is still clearly visible, over at least one order of magnitude, even in the presence of the steepest trap (with σ trap, f = 100 in figure 5).
One may ask about the role of interactions between weakly interacting atoms of a Bose condensate. Those effectively scale in a mean-field approximation with the number of atoms, which may be controlled and possibly be reduced [14]. We modeled the evolution of an initially well-confined condensate, following the above protocol, using a 1D Gross-Pitaevskii equation. In this approach, the interactions are taken into account by a nonlinear density-dependent term in the Hamiltonian [14,34,35]. A weak repulsive nonlinearity may actually stabilize the evolution and-to some extent-reverse the effect of a weak confinement during the  (21), for a direct comparison with figure 3 and the power law of equation (7), with (α − 1)/2 = 1.5 (thin solid line). relaxation (cf [37] for a similar effect). This can be seen in our figure 6, where the prefactor of the nonlinear term is denoted by g in our dimensionless units 8 .

Conclusions
In summary, we propose a setup for which one may observe a power-law decay of the asymptotic survival probability P S (t) in a controlled manner. Such an anomalous decay is readily realized by preparing a quantum wave packet with a sufficiently broad energy range around the ground state energy E = 0 in the potential wells given by equations (1) and (19). The slow algebraic 8 Our dimensionless potentials are given in units of energy m −2 . The three-dimensional atom-atom interaction strength can be reduced to an effective 1D parameter, provided a strong transverse/radial confinement is experimentally achieved [14]. The latter parameter g 1D = 2hω rad a S N a , with the scattering length a S , the radial confinement frequency ω rad and the number of atoms in the Bose condensate N a , can be expressed without dimensions using, e.g., just the single scale given by the radial confinement: E 0 =hω rad and x 0 = √ 2/hω rad M, where M is the single atom mass. This gives g = g 1D /(E 0 x 0 ) = 2N a a S /x 0 . Here the wave function is normalized to 1 as in our single-particle computations otherwise used in the paper. relaxation toward this equilibrium state arises from the spectral properties and the population of eigenmodes of the system. The latter are encoded in the coefficient C 2 defined in equation (5). Hence, the predicted behavior is of purely quantum origin (not induced, however, by dynamical or Anderson localization [4,5] as discussed in [38], nor by semiclassical arguments [39]). The exceptional control of state-of-the-art experiments with ultracold atoms [15,[18][19][20][21]36] offers the possibility to observe our predictions following, e.g., our protocol for obtaining the data of figures 5 and 6.
with α = √ 1 + 8V 1 − 1. By imposing the continuity in x = ±L, one obtains By imposing that the ground state is normalized to unity, one obtains The continuity of the first derivative in x = ±L is ensured by equation (2). The odd eigenfunctions for the continuum part of the spectrum are where ν = (α + 1)/2. By imposing the continuity in x = ±L, one obtains By imposing the normalization condition dx ψ E (x)ψ E (x) = δ(E − E ), one obtains a E = 1 √ 2 cos( ν + (E)), b E = 1 √ 2 sin( ν + (E)). (A.8) By imposing the continuity of the first derivative in x = ±L, one obtains