From the area under the Bessel excursion to anomalous diffusion of cold atoms

Levy flights are random walks in which the probability distribution of the step sizes is fat-tailed. Levy spatial diffusion has been observed for a collection of ultra-cold Rb atoms and single Mg+ ions in an optical lattice. Using the semiclassical theory of Sisyphus cooling, we treat the problem as a coupled Levy walk, with correlations between the length and duration of the excursions. The problem is related to the area under Bessel excursions, overdamped Langevin motions that start and end at the origin, constrained to remain positive, in the presence of an external logarithmic potential. In the limit of a weak potential, the Airy distribution describing the areal distribution of the Brownian excursion is found. Three distinct phases of the dynamics are studied: normal diffusion, Levy diffusion and, below a certain critical depth of the optical potential, x~ t^{3/2} scaling. The focus of the paper is the analytical calculation of the joint probability density function from a newly developed theory of the area under the Bessel excursion. The latter describes the spatiotemporal correlations in the problem and is the microscopic input needed to characterize the spatial diffusion of the atomic cloud. A modified Montroll-Weiss (MW) equation for the density is obtained, which depends on the statistics of velocity excursions and meanders. The meander, a random walk in velocity space which starts at the origin and does not cross it, describes the last jump event in the sequence. In the anomalous phases, the statistics of meanders and excursions are essential for the calculation of the mean square displacement, showing that our correction to the MW equation is crucial, and points to the sensitivity of the transport on a single jump event. Our work provides relations between the statistics of velocity excursions and meanders and that of the diffusivity.


INTRODUCTION
The velocity, v(t), of a particle interacting with a heat bath exhibits stochastic behavior which in many cases is difficult to evaluate. The position of the particle, assumed to start at the origin at time t = 0, is the time integral over the fluctuating velocity x(t) = t 0 v(t )dt and demands a probabilistic approach to determine its statistical properties. Luckily the central limit theorem makes it possible, for many processes, to predict a Gaussian shape for the diffusing packet. Then the diffusion constant K 2 characterizes the normal motion through its mean square displacement x 2 = 2K 2 t. The remaining goal, for a given process or model, is to compute K 2 and other transport coefficients. In normal cases this can be done, at least in principle, via the Green-Kubo formalism, namely by the calculation of the stationary velocity correlation function, which gives the diffusivity An alternative approach is investigated in this work and is based on the concept of excursions. We assume that the random process v(t) is recurrent and thus the velocity crosses the zero point, v = 0, many times in the observation window (0, t). We divide the path x(t) into a sum of increments  (1) Here {t 1 , t 2 · · ·} are the points in time of the velocity zero-crossings, v(t i ) = 0. In the interval (t i , t i+1 ) the velocity is either strictly positive or negative. The velocity in each interval is thus a stochastic process which starts and ends on the origin without crossing it in between. Such a random curve is called an excursion. The random spatial increment χ i = ti+1 ti v(t )dt is the area under the excursion. The position of the particle, according to Eq.
(1), is the sum of the random increments, namely a sum of the signed areas under the velocity excursions, each having a random duration. The goal of this paper is to relate the statistics of the areas under these velocity excursions and the corresponding random time intervals between zero crossings, to the problem of spatial diffusion. This connection is easy to find in the case that the increments χ i and the duration of the excursions τ i = t i+1 −t i are mutually uncorrelated, independent and identically distributed random variables. Over a long measurement time, the number of excursions is t/ τ where τ is the average time for an excursion. Then according to Eq. (1) the mean squared displacement is x 2 = χ 2 t/ τ , and hence we have This equation is reminiscent of the famous Einstein formula, and shows that diffusion is related to the statistics of excursions. The original work of Einstein, discussed in many textbooks, is explicitly based on an underlying random walk picture, and so does not involve the area under random velocity excursions, nor zero crossings in velocity space. We will arrive at the simple equation (2) only at the end of our work, in Sec. , while here it merely serves as an appetizer to motivate the consideration of velocity excursions in some detail, and to suggest the usefulness of developing tools for the calculation of χ 2 and τ (here χ = 0, by symmetry). Obviously Eq. (2) is based on the assumption that the variance χ 2 is finite, and as mentioned, that correlations are not important. The major effort of this paper is directed to considering a more challenging case, namely a physically relevant stochastic process where the variance χ 2 diverges, and more importantly the process exhibits correlations between χ and τ . The particular system we investigate is a model for diffusion of atoms in optical lattices, following the experiments [1][2][3] and the semi-classical theory of Sisyphus cooling [4][5][6]. One type of excursion which has been thoroughly investigated is the Brownian excursion [7][8][9]. A Brownian excursion is a conditioned one dimensional Brownian motion v(t) over the time interval 0 < t < τ . The motion starts at v(0) = (eventually the → 0 limit is taken) and ends at v(τ ) = and is constrained not to cross the origin v = 0 in the observation time (0, τ ). The area under this curve is a random variable, whose statistical properties have been investigated by mathematicians [10][11][12][13]. More recently, this problem was treated with a path integral approach to describe statistics of fluctuating interfaces [7,8,14,15] and related areas until the first passage time are used to describe universality of sandpile models [16]. Here, we generalize the Brownian excursion to a process v(t) described by a Langevin equation, with an asymptotically logarithmic potential. We show how the random area under this Langevin excursion determines the dynamics of cold atoms. We believe that Langevin excursions, or more generally random excursions, are useful tools in many areas of statistical physics, hence their investigation beyond the well studied Brownian excursion, is worthwhile.
If v(t) is a Brownian motion, or as we will proceed to show for cold atoms in shallow lattices, the set of points where v(t) = 0 is non-trivial, as it contains no isolated points and no interval. This is the case since we are dealing with a continuous path with power law statistics of the crossing times (see details below). Mathematician have investigated the statistics of level crossing, e.g. zero crossing, of continuous paths in great detail. The concept of local time, introduced by P. Lévy in the context of Brownian motion and Ito's excursion theory for continuous paths, are pillars in this field, see [17,18] and references within. Here we use a heuristic approach, with a modification of the original process v(t), by introducing a cutoff ± : the starting point of the particle after each zero hitting, which is taken to zero at the end. This makes it possible to use renewal theory, and continuous time random walks, which are tools well studied in physics literature of discrete processes. In this sense we differ from rigorous mathematical approaches. Thus we avoid the problem of continuous paths, for example the infinite number of zero crossings, by replacing the original path with an modified path for which the number of zero crossings is finite (when is finite). We show that physical quantities characterizing the entire process have a finite → 0 limit. For example in Eq. (2), for Langevin dynamics of v(t), both τ and χ 2 approach zero as → 0, their ratio K 2 approaches a finite limit.
One might wonder why we wish to use excursions and their peculiar properties to evaluate the spreading of atoms or more generally other transport systems. The answer is that it turns into a useful strategy when the friction forces are non-linear. In particular we will investigate laser cooling, where within the semi-classical theory, the dimensionless friction force is [5] (see details below) This friction force, induced by the laser fields, is linear for small velocities F (v) ∼ −v similar to the Stokes friction law for a massive Brownian particle in water at room temperature. However, unlike such friction, which increases in magnitude with velocity, here for large v, F (v) ∼ −1/v → 0. Asymptotically, then, the system is frictionless. This implies that fast particles have to remain fast for a long time, which in turn induces heavytailed populations of fast particles [19]. In this case, as we show later, the standard picture of diffusion breaks down. More specifically, the problem of diffusion of cold atoms under Sisyphus cooling was partially treated by Marksteiner, et al. [6]. While clearly indicating the anomalous nature of the diffusion process, the main tool used was the evaluation of the stationary velocity correlation function of the process, followed by the use of the Green-Kubo formalism for the evaluation of K 2 . They showed that for a certain critical value of the depth of the optical lattice the value of K 2 diverges (see also [20]). Katori et al. [1] measured the mean square displacement x 2 ∼ t 2ξ and recorded the onset of super-diffusion ξ > 1/2 beyond a critical depth of the optical lattice [1]. Wickenbrock et al. [2] in the context of driven optical lattice experiments demonstrated with Monte Carlo simulations an upper limit on the spreading of the atoms x 2 ≤ const × t 3 . Sagi et al. [3] showed that a packet of spreading Rb atoms can be fitted with a Lévy distribution instead of a Gaussian distribution found for normal diffusion. These findings clearly indicate the breakdown of the usual strategy of treating normal diffusion and hence promoted further theoretical investigations. The velocity correlation function is not stationary and hence the Green-Kubo formalism must be replaced [22]. The moments exhibit multifractal behavior [23], there is an enhanced sensitivity to the initial preparation of the system [24,25], momentum fluctuations are described by an infinite covariant density [26] and in certain parameter regimes the Lévy central limit theorem applies instead of the standard Gaussian version [6,21]. In this regime of shallow optical lattices the power of the analysis of the area under Langevin excursions becomes essential as we will show here. We note that the relation of laser cooling with Lévy statistics is not limited to the case of Sisyphus cooling considered in this manuscript. Sub-recoil laser cooling, a setup different from ours, also leads to fundamental relations between statistical physics of rare events, and laser-atom physics [27,28]. For a recent minireview on the departure from Boltzmann-Gibbs statistical mechanics for cold atoms in optical lattices, see [29].

Scope and organization of paper
The current work significantly extends the investigation of the properties of the spatial distribution of atoms in Sisyphus cooling begun in Ref. [21]. There we uncovered three phases of the motion which are controlled by the depth U 0 of the optical lattice: a Gaussian phase x ∼ t 1/2 , a Lévy phase x ∼ t ξ and 1/2 < ξ < 3/2, and a Richardson phase x ∼ t 3/2 (see details below). Within the intermediate phase, the density of particles, in the central region of the packet, is described by a symmetric Lévy distribution, similar to the fitting procedure used in a recent Weizmann Institute experiment [3]. However, this cannot be the whole story. As is well known the variance of the Lévy distribution diverges, which implies that x 2 = ∞, which is unphysical. Indeed, as mentioned, Katori et al. experimentally determine a finite mean square displacement x 2 ∼ t 2ξ (see also [2]). We showed related numerical evidence that the Lévy distribution is cut off at distances of the order x ∼ t 3/2 . This breakdown of Lévy statistics arises due to the importance of correlations between jump lengths χ and jump duration τ , which are neglected in the derivation of the Lévy distribution. Our purpose here is to investigate these correlations in detail. As we proceed to show, the correlations are given by the statistical properties of the Langevin excursions discussed above, for all except the last interval. Because the velocity is not constrained to be zero at the measurement time t, this last interval is not described by an excursion but rather by a Langevin meander, where the random walk v(t) begins at the origin and does not return for the entire duration of the walk. The properties of the excursions and meanders enter in a modified Montroll-Weiss [30] equation for the Fourier-Laplace transform of the density P (x, t) which we derive. We then use this equation to calculate a quantity which is sensitive to the correlations, namely the mean square displacement. The mean square displacement exhibits anomalous diffusion and is sensitive to the last jump event, i.e., to the statistics of the meander. Thus, our treatment modifies both the celebrated Montroll-Weiss equation, to include the last jump in the sequence (i.e., the meander) and the existing theory of areas under Brownian meanders and excursions to include the dissipative friction force, which is responsible to the cooling of the atoms. Our analysis illuminates the rich physical behavior and provides the needed set of mathematical tools beyond the decoupling approximation used to obtain the Lévy distribution in our previous work. For completeness, we present a detailed decoupled and coupled analysis, the latter being the main focus of the current work.
The paper is organized as follows. We start with a brief survey of the semi-classical theory of Sisyphus cooling [5,6], and show the connection of the dynamics to Lévy walks following Marksteiner et al. [6]. The importance of the correlations between τ and χ is emphasized, a theme which as mentioned has not received its deserved attention. In Sec. , a simple scaling theory is presented which yields the exponents describing the dynamics of the atomic packet. The main calculation of the distribution of the area under the Bessel excursion is found in Sec. , the calculation of the area under the Bessel meander is given in an Appendix. A new coupled continuous time random walk theory in Sec. provides the connection between the statistics of excursions and meanders, and the evolution of the density profile. Asymptotic behaviors of the Fourier-Laplace transform of the joint probability density function (PDF) of jump lengths and waiting times are investigated in Sec. . These in turn give us the asymptotic behaviors of the atomic density packet, the mean square displacement, and the different phases of the dynamics which are investigated in Secs.
-. Derivation of the distribution of the time interval straddling time t for the Bessel process is carried out in Appendix F, which allows as to connect between our heuristic renewal approach and more rigorous treatments [17,31] and further discuss the nontrivial fractal set of zero crossings.

SEMI-CLASSICAL DESCRIPTION OF COLD ATOMS-LANGEVIN DYNAMICS
We briefly present the semi-classical picture for the dynamics of the atoms. The trajectory of a single particle with mass m is x(t) = t 0 p(t)dt/m where p(t) is its momentum. Within the standard picture [5,6,33] of Sisyphus cooling, two competing mechanisms describe the dynamics. The cooling force F (p) = −αp/[1 + (p/p c ) 2 ] acts to restore the momentum to the minimum energy state p = 0. Momentum diffusion is governed by a diffusion coefficient which is momentum dependent, D(p) = . The latter describes momentum fluctuations which lead to heating due to random emission events which stochastically jolt the atom. We use dimensionless units, time t → tα, momentum p → p/p c , distance x → xmα/p c , and introduce the dimensionless momentum diffusion constant D = D 1 /(p c ) 2 α. For simplicity, we set D 2 = 0 since it does not modify the asymptotic |p| → ∞ behavior of the diffusive heating term, nor that of the force and therefore does not modify our main conclusions.
The Langevin equations dp dt describe the dynamics in phase space.
Here the noise term is Gaussian, has zero mean and is white, The force F (p) is a very peculiar non-linear friction force. We note that friction forces which decrease with increasing velocity or momentum like 1/p are found also for some nanoscale devices, e.g., an atomic tip on a surface [32]. One goal of the paper is to find the spatial distribution of particles governed by Eq. (4). The stochastic Eq. (4) gives the trajectories of the standard Kramers picture for the semi-classical dynamics in an optical lattice which in turn was derived from microscopic considerations [5,6]. Denoting the joint PDF of (x, p) by W (x, p, t), the Kramers equation reads From the semiclassical treatment of the interaction of the atoms with the counter-propagating laser beams, we have where U 0 is the depth of the optical potential, E R the recoil energy and the dimensionless parameter c [34] depends on the atomic transition involved [5,6,35]. U 0 is a control parameter; hence different values of D are attainable in experiment, and exploration of different phases of the dynamics are within reach [1][2][3]19]. Eq. (7) is rather intuitive since deep optical lattices, i.e. large U 0 , implies small D while large recoil energy leads to a correspondingly large value of D.
The behavior of the distribution for momentum only (when x is integrated out from the Kramer's Eq., yielding the Fokker-Planck Eq.) is much simpler, and has been presented in previous work [26]. The equilibrium properties are governed by the effective potential in momentum space, which for large p is logarithmic, V (p) ∼ ln(p). This large p behavior of V (p) is responsible for several unusual equilibrium and non-equilibrium properties of the momentum distribution [1,19,24,26,36]. The equilibrium momentum distribution function is given by [19] W eq (p) Here (10) is the normalizing partition function. Eq. (9) is Student's t distribution, also sometimes called a Tsallis distribution [19]. Actually the problem is related coincidentally to Tsallis statistics, since as mentioned the equilibrium PDF is proportional to a Boltzmann-like factor, W eq (p) ∝ exp[−V (p)/D], so D acts like a temperature. More importantly, the power-law tail of the equilibrium PDF, for sufficiently large D, implies a large population of fast particles, which in turn will spread in space faster than what one would expect using naive Gaussian central limit theorem arguments. For example if 1/D < 3 the ensemble averaged kinetic energy in equilibrium diverges, since p 2 eq = ∞ while when 1/D < 1 the partition function diverges and a steady-state equilibrium is never reached. A dramatic increase of the energy of atoms when the optical lattice parameter U 0 approaches a critical value was found experimentally in [1] and a power-law momentum distribution was measured in [19]. Of course, the kinetic energy of a physical system cannot be infinite, and the momentum distribution must be treated as a time dependent object within the infinite covariant density approach [26]. While much is known about the momentum distribution, both experimentally and theoretically, the experiments [1-3] demand a theory for the spatial spreading. We note that diffusion in logarithmic potentials has a vast number of applications [24][25][26][37][38][39][40][41] and refs. therein, e.g. Manning condensation [42], and unusual mathematical properties, including ergodicity breaking [36]. In general, logarithmic potentials play a special role in statistical physics [43][44][45].

MAPPING THE PROBLEM TO A LÉVY WALK PROCESS
In principle one may attempt to directly solve the Kramers equation (6) to find the joint PDF of the random variables (x, p) at a given time t. In Fig. 1 we plot a histogram of the phase space obtained numerically. We see a complicated structure: roughly along the diagonal, clear correlations between x and p are visible; on the other hand, along the x and p axis, decoupling between momentum and position is evident, together with broad (i.e. non-Gaussian) distributions of x and p. At least to the naked eye no simple scaling structure is found in x−p space, and hence we shall turn to different, microscopic, variables which do exhibit simple scaling. This leads to an analysis centered on the mapping of the Langevin dynamics to a Lévy walk scheme [6,46] and the statistics of areas under random excursions.
Starting at the origin, p = 0, the particle along its stochastic path in momentum space crosses p = 0 many times (see Fig. 2). In other words the random walk in momentum space is recurrent; this being the case even when D → ∞ since then the process in momentum space is one of pure diffusion (i.e. the force is negligible) and from Polya's theorem we know that such one dimensional walks are recurrent. The cooling force being attractive clearly maintains this property.
Let τ > 0 be the random time between one crossing event and the next and let −∞ < χ < ∞ be the random displacement for the corresponding τ . As shown schematically in Fig. 2 the process starting with zero momentum is defined by a sequence of jump durations {τ 1 , τ 2 , · · ·} with corresponding displacements {χ 1 , χ 2 , · · ·}. These random waiting times and displacements generate a Lévy walk [6], as we explain below. Here p(t )dt etc. Let the points on the time axis {t 1 , t 2 , · · ·} denote times t n > 0 where the particle crossed the origin of momentum p = 0 (see Fig. 2). These times are related to the waiting times: dt is the area under the random momentum curve constrained in such a way that p(t ) in the time interval (t k−1 , t k ) does not cross the origin, while it started and ended there. The total displacement, namely the position of the particle at time t, is a sum of the individual displacements x = n i=1 χ i +χ * . Here n is the random number of crossings of zero momentum in the time interval (0, t) and χ * is the displacement made in the last interval (t n , t). By definition, in the time interval (t n , t) no zero crossing was recorded. The time τ * = t − t n is sometimes called the backward recurrence time [47]. The measurement time is clearly t = n i=1 τ i + τ * . In standard transport problems the effect of the last displacement χ * on the position of particle x is negligible, and similarly one usually assumes t t n when t is large. However, for anomalous diffusion of cold atoms, where the distributions of displacements and jump durations are wide, these last events cannot be ignored.
One goal is to find the long time behavior of P (x, t), the normalized PDF of the spatial position of a particle that started at the origin x = 0, p = 0 at t = 0. It is physically clear that in the long time limit our results will not be changed if instead we consider narrow initial conditions, for example Gaussian PDFs of initial position and momentum. Initial conditions with power-law tails will likely lead to very different behaviors [24,25]. Once we find P (x, t) we have the semi-classical approximation for the spatial density of particles. The latter can be compared with the Weizmann Sisyphus cooling experi-ments [3] provided that collisions among the atoms are negligible.
The Lévy walk process under investigation is a renewal process [47] in the sense that once the particle crosses the momentum origin the process is renewed (since the Langevin dynamics is Markovian). This is crucial in our treatment, and it implies that the waiting times τ i are statistically independent identically distributed random variables as are the χ i . However, as we soon discuss, the pairs {τ i , χ i } are correlated.
Since the underlying Langevin dynamics is continuous, we need a refined definition of the Lévy walk process. Both the τ i 's and the χ i 's are infinitesimal; however the number of renewal events, n, diverges for any finite measurement time t, in such a way that the total displacement x is finite. In this sense the Lévy walk process under investigation is different from previous works where the number of renewals, for finite measurement time is finite. One way to treat the problem is to discretize the dynamics, as is necessary in any event to perform a computer simulation, and then χ i and τ i are of course finite. In our analytical treatment, following Marksteiner, et al. [6], we consider the first passage time problem for a particle starting with momentum p i and reaching p f < p i for the first time at τ . We take p f = 0 and eventually take p i = → 0. The Lévy walk scheme is hence summarized with the following steps: 1. Choose with probability 1/2 either +p i or −p i .
2. Follow the Langevin dynamics until the particle reaches p f = 0.
3. Record the random displacement χ and random duration τ during this excursion.

Go to 1.
This loop is terminated at time t, the final displacement χ * calculated, and as mentioned the total displacement is x = n i=1 χ i + χ * . In the first step we have probability 1/2 to start with either +p i or −p i since the cooling force is antisymmetric and so vanishes at p = 0. The advantage of presenting the problem as a set of recurrent random walks through the χ's and τ 's instead of the direct Langevin picture stems from the fact that we can treat analytically the former Lévy walk picture.
We denote the joint PDF of the pair {χ, τ } of a single excursion by ψ(χ, τ ). The theoretical development starts from the analysis of ψ(χ, τ ) and then from this we use the Lévy walk scheme to relate this single excursion information to the properties of the entire walk, and in particular, P (x, t) for large t.
There exists a strong correlation between the excursion duration τ and the displacement length χ which is encoded in ψ(χ, τ ). Let the PDFs of χ and τ be denoted q(χ), g(τ ), respectively. Weak correlations would imply the decoupled scheme ψ(χ, τ ) q(χ)g(τ ) which We plot the flight distance |χ| versus the jump duration τ to demonstrate the strong correlations between these two random variables. Here D = 2/3 and the red line has a slope 3/2 reflecting the χ ∼ τ 3/2 scaling discussed in the text.
is far easier to analyze; however as we shall see this decoupling is not generally true and leads to wrong results for at least some observables of interest. Properties of q(χ) and g(τ ) were investigated in [6,35] and are also studied below. The problem becomes more interesting and challenging due to these correlations. As we show, these are responsible for the finiteness of the moments of P (x, t), in particular the mean square displacement x 2 , and for the existence of a rapidly decaying tail of P (x, t). This in turn is related to the Lévy flight versus Lévy walk dilemma [46], to multifractality [23], and to the physical meaning of the fractional diffusion equation [48] used as a fitting tool for the Weizmann experiment [3] (see Eq. (105) and discussion there). As we show below, beyond a critical value of D = 1 the correlations can never be ignored and govern the behavior of the entire packet, not only the large x tails of P (x, t). Physically, the correlations are obvious, since long durations of flight τ involve large momentum p, which in turn induce large displacements χ. As an example of these correlations, we plot in Fig. 3 the displacement |χ| versus the corresponding τ obtained from computer simulation. The figure clearly demonstrates the correlations and it also shows a |χ| ∝ τ 3/2 scaling which we now turn to investigate. Notice should be paid to how much simpler the χ−τ distribution is, compared to the x−p distribution in Fig. 1. Simulations presented in Fig. 3 were performed on a discrete lattice in momentum space, starting on the first lattice point, see Appendix for details.

SCALING THEORY-RELATION BETWEEN EXPONENTS
As shown by Marksteiner, et al. [6], and Lutz [35] the PDFs of the excursion durations and displacements satisfy the asymptotic laws In Appendices A and B we study these PDFs using backward Fokker-Planck equations. In particular we find the amplitudes g * and q * , and the relevant PDFs moments.
What is most crucial are the exponents β and γ. When D → ∞ we get γ = β = 0 and Eq. (11) gives familiar limits. In the "high temperature" limit of large D, the cooling force is negligible and then the Langevin equation reduces to Brownian motion in momentum space. Then g(τ ) ∝ τ −3/2 , which is the well-known asymptotic behavior of the PDF of first passage times for unbounded one dimensional Brownian motion [49,50]. Less well known is lim D→∞ q(χ) ∝ |χ| −4/3 which describes the distribution of the area under a Brownian motion until its first passage time (see [51,52] who give this PDF explicitly). Notice that the power-law behavior Eq. (11) yields a diverging second moment of the displacement χ for D > 1/5 which in turn gives rise to anomalous statistics for x. The correlations between χ and τ are now related to the asymptotic behaviors of their PDF's, Eq. (11). We rewrite the joint PDF ψ(χ, τ ) = g(τ )p(χ|τ ) (13) where p(χ|τ ) is the conditional PDF to find jump length χ for a given jump duration τ . We introduce a scaling ansatz which is expected to be valid at large τ : Since Changing variables to z = χ/τ η , we get Comparing with Eq. (11) we get a simple equation for the unknown exponent η which is 1 + 1/(2η) + γ/η = 4/3 + β, so using Eq. (12) we find η = 3/2. This is precisely the scaling behavior χ ∼ τ 3/2 we observe in our simulation, Fig. 3. Hence the natural scaling solution of the problem is In hindsight this result is related to Brownian scaling. For a particle free of the cooling force, its momentum will exhibit Brownian scaling p ∼ (Dτ ) 1/2 and hence the excursions which are integrals of the velocity scale as χ ∼ √ Dτ 3/2 . The crucial point is that this simple Brownian scaling is maintained even in the presence of the cooling force for all D. This is due to the marginal nature of the weak 1/p friction force for large p which leaves the Brownian scaling intact, but changes the scaling function. While the 3/2 scaling is simple and D independent, the shape of B(·) is sensitive to this control parameter. In the next section we investigate B(·) and for this we must go beyond simple scaling arguments.

AREA UNDER THE BESSEL EXCURSION
A natural generalization of the Brownian excursion is a Langevin excursion. Such a stochastic curve is the path p(t ), given by the Langevin equation, in the time interval 0 ≤ t ≤ τ , such that it starts and ends at p i = p f = , but is constrained to remain positive in between. Here, p i = p(0) is the initial and p f = p(τ ) is the final location in momentum space. For our application, the path is considered in the limit → 0. Since the path never crosses the origin, the area under such a curve is χ = τ 0 p(t )dt , and hence the PDF of χ for fixed τ yields the sought after conditional PDF, p(χ|τ ). Obviously, χ is the integral over the constrained path p(t ); hence by definition it is the area under the excursion. The meander will describe the last jump χ * since at the measurement time the particles velocity is generally non-zero. For now we will discuss only excursions, and find p(χ|τ ), and return to the meander later.
Here we focus our attention on a specific excursion we call the Bessel excursion, corresponding to the case F (p) = −1/p: so that the effective potential is the non-regularized logarithm, V (p) = ln(p) and p > 0. Since the scaling approach is valid for long times, where excursions are long, the typical momentum p is large and the details of the force field close to the origin are negligible for the purpose of the calculation of the scaling function B(·). We will check this assumption with numerical simulations, which of course use a regularized form of the force law. Some sample paths of Bessel excursion are presented in Fig.  4. The name Bessel excursion stems from the fact that Langevin dynamics in the non-regularized friction force The random excursions are constrained Langevin paths, with a −1/p force, that do not cross the origin in the observation time t while start and end on → 0. For simulations we used the regularized force field, which only alters the dynamics when p 1, and is negligible in the the long time limit (see Appendix D). When D → ∞ we get a Brownian excursion. We see that as D is decreased the excursions are further pushed from the origin p = 0, since small D implies effectively large forces, hence to avoid the zero crossing the particles must drift further away from origin. Thus the attractive force, repels particles, which at first might sound counter intuitive and seems to be an unexplored property of excursions. field 1/p corresponds to a well-known stochastic process called the Bessel process [41,53]. More information on the regularized and non-regularized processes is given in Appendixes A,B.
Let G τ (χ, p|p i ) be the joint PDF of the random variables χ and p. Since Later we will take p to be the final momentum p f , which similarly to p i will be set to the value → 0 (for the sake of notational brevity we omit the subscript in p f ). The calculation of p(χ|τ ) follows three steps. For the Brownian case F (p) = 0, this method was successfully applied by Majumdar and Comtet [8].
The reason why we multiply s with D in the Laplace transform will become clear soon. Since χ is a functional of the path p(t ) we will use the Feynman-Kac (FK) formalism to find G τ (s, p|p i ) (see details below). The constraint that the path p(t ) is always positive enters as an absorbing boundary condition at p = 0 [41].
(ii) The second step is to consider p i = p = → 0 and obtain the Laplace transform The denominator in the above equation ensures the normalization, since we must havep(s|τ )| s=0 = 1.
We now implement these steps to solve the problem. The FK formalism [52] treats functionals of Brownian motion. Here we use a modified version of the FK equation to treat over-damped Langevin paths [54].
dt be a functional of the Langevin path and assume U (·) > 0 so that we are treating positive functionals. Here G τ (A, p|p i ) is the joint PDF of A and p and G τ (s, p|p i ) is the corresponding Laplace A → sD transform. The generalized FK equation reads is the celebrated FK equation which is an imaginary time Schrödinger equation and −sDU (p) is the potential of the corresponding quantum problem. The constraint on positive excursions, namely p > 0, gives the boundary condition G τ (s, 0|p i ) = 0. In the quantum language this is an infinite potential barrier for p < 0. This formalism can be used in principle to obtain the area under Langevin excursions for all forms of F (p). For the Bessel excursion under investigation here, the functional U [p(t)] = p(t) is linear since χ = τ 0 p(t )dt . Quantum mechanically, this gives a linear potential and hence the connection to the Airy function found for F (p) = 0 in [7,8]. With the force field F (p) = −1/p we haveL fp = D(∂ p ) 2 + ∂ p p −1 and hence we find, using Eqs. (18,19,21) (the former gives the sD term) The solution of this equation is found using the separation ansatz which yields the time independent equation We now switch to the more familiar one dimensional Schrödinger equation via the similarity transformation This Schrödinger equation has a binding potential, which yields discrete eigenvalues, with an effective repulsive potential with a p −2 divergence for p → 0 and a binding linear potential for large p provided s = 0. Eq. (25) describes a three dimensional non-relativistic quantum particle in a linear potential [55,56], the p −2 part corresponding to an angular momentum term. As usual φ k (p) yields a complete orthonormal basis ∞ 0 φ k (p)φ m (p)dp = δ km and the formal solution of the problem is We can scale out the Laplace variable s from the time independent problem, defining φ k (p) = s 1/6 f k (s 1/3 p) and  It follows from Eq. (27) and the absorbing boundary condition that withd k a k-dependent coefficient. Thed k will soon be seen to be important and they are evaluated from solutions of Eq. (27) with the normalization condition ∞ 0 [f k (p)] 2 dp = 1. For the initial and final conditions under investigation, p i = p = , we have where the exponents ν and α will turn out to be useful Note that classification of boundaries for a nonregularized Bessel process was carried out in [41] and discussed here in Sec. . For an absorbing boundary condition both the sign of the probability current and usual condition of the vanishing of the probability on the absorbing point must be taken into consideration [41]. According to Eq. (20), we need G τ (s = 0, | ), in order to normalize the solution. This s = 0 propagator can be found exactly. The eigenvalue problem now reads The superscript 0 indicates the s = 0 case. Since s = 0 the linear field in Eq. (25) is now absent so the "particle" is not bounded and hence one finds a continuous spectrum E 0 k = k 2 . The wave functions consistent with the boundary condition are where J α (·) is the Bessel function of the first kind. The second solution with J −α (kp) is unphysical due to the boundary condition [41]. The normalization condition is where L → ∞ is the "box" size which eventually will drop out of the calculation. Solving the integral in Eq. (33) gives B k = πk/L. For initial and final conditions p f = p i = the s = 0 propagator then reads Using the small z expansion of the Bessel function Notice that replacing J α (·) with J −α (·) leads to G τ (s = 0, | ) ∼ 1−2α which diverges since α > 1/2, as mentioned this unphysical solution is rejected due to the absorbing boundary conditions [41]. The usual density of states calculation follows the quantization rule of a particle in a box extending between 0 and L which gives kL = nπ with integer n. Hence k · · · → ∞ 0 dk L π · · · and one finds, after a change of variable to x = k 2 , Inserting Eqs. (29,36) in Eq. (20) we find our first main resultp In the limit D → ∞ we find |d k | → 1 and the λ k 's are the energy eigenvalues of the Airy equation, related to the zeros of the Airy function, and, up to a rescaling of s, we get the result obtained by Darling [10] and Louchard [11] for the area under the Brownian excursion, namely the Laplace transform of the Airy distribution [8]. We can show thatp(s = 0|τ ) = 1 as it should. It is easy to tabulate the eigenvalues λ k and the coefficientsd k using Eq. (27) and standard numerically exact techniques. In Table I we tabulate the first few slopesd k and eigenvalues for D = 2/5. For large k, i.e. large energy, the 1/p 2 part of the potential is irrelevant and the eigenvalues λ k converge to the eigenvalues of the Airy problem considered previously. Similarly, we expect lim k→∞ |d k | = 1 for any finite D, since in the Airy problem limit, i.e., D → ∞ case [7,8], |d k | is unity for all k.
To complete the calculation we need to perform an inverse Laplace transform, namely invert from sD back to χ > 0 (and multiply the PDF by 1/2 if interested in both positive and negative excursions). In Eq. (37) we have terms with the structure Each term on the right hand side can be inverse transformed separately since the inverse Laplace transform of (Ds) γ is χ −γ−1 /Γ(−γ). After term by term transformation we sum the infinite series (summation over n) using Maple. Thus we arrive at our first destination, the conditional PDF for χ > 0 is found in terms of generalized hypergeometric functions where the summation is over the eigenvalues. In Fig. 5 we plot the solution for D = 2/5 and D → ∞ corresponding the the Brownian case. Notice the χ ∼ √ Dτ 3/2 scaling which proves the scaling hypothesis Eq. (17). The sum in Eq. (39) converges quickly as long as χ is not too large, and so we can use it to construct a plot of p(χ|τ ). For example, for Fig. 5 only k = 0, ...4 is needed to obtain excellent convergence. We stress that Eq. (39) gives an explicit representation of the scaling function for 0 < χ < ∞, and hence by symmetry for all χ. The scaling variable is and the scaling function B(v 3/2 ) in Eq. (17) can be read directly off of Eq. (39).

The Bessel Meander
As noted, the last zero crossing of the momentum process p(t) takes place at a random time t n , and hence at time t the particle is unlikely to be on the origin of momentum. Since the particle, by definition, did not cross the origin in (t n , t) its momentum remains either positive or negative in the backward recurrence time interval τ * = t − t n (with equal probability). This means that in this time interval the motion is described by a meander, not an excursion. The area under the Brownian meander was investigated previously [8,9,57]. For our purposes we need to investigate the Bessel meander. This is a Langevin path described by Eq. (18) constrained to remain positive, which starts at the origin but is free to end with p > 0. More specifically, the diffusive scaling of χ * ∼ (τ * ) 3/2 still holds. Then as for the pairs (χ, τ ), we have the conditional PDF Here the subscript M stands for a meander. The scaling function B M (.) is different from B(.) though clearly both are symmetric, with mean equal to 0 (since positive and negative meanders and excursions are equally probable). The calculation of B M (.) runs parallel to that for the excursion and is presented in Appendix E. Soon we will demonstrate the importance of the meander for specific observables of interest. Having explicit information on B(.) and g(τ ), Eq. (11) (see Appendix A for details) and the scaling form of B M (.) (see Appendix E), we can now investigate the packet P (x, t).

MONTROLL-WEISS EQUATION FOR FOURIER-LAPLACE TRANSFORM OF P (x, t)
We now use tools developed in the random walk community [58][59][60][61] to relate the joint PDF of a single excursion Eqs. (13,14), to the probability density P (x, t) for the entire walk. We find a modified Montroll-Weiss [30,48,62] type of equation for the Fourier-Laplace transform of P (x, t) in terms of the Fourier-Laplace transform of ψ(χ, τ ) which will be denotedψ(k, u). One modification is that we include here the correct treatment of the last jump. Usually, the continuous time random walk (CTRW) model [48], has as an input a single joint PDF of jump lengths and times, while in our case we have essentially two such functions describing the excursions (i.e. B(.)) and the meander (i.e. B M (.)). Generally, Montroll-Weiss equations are the starting point for derivation of fractional diffusion equations and the asymptotic behaviors of the underlying random walks [48]. The original work of Montroll and Weiss [30] assumed there were no correlations between step size χ and and the waiting time τ , corresponding to a situation called a decoupled CTRW. The diffusion of atoms in optical lattices corresponds to a coupled spatial-temporal random walk theory first considered by Scher and Lax [63] (see [64][65][66] for recent developments). Define η s (x, t)dtdx as the probability that the particle crossed the momentum state p = 0 for the sth time in the time interval (t, t + dt) and that the particle's position was in the interval (x, x + dx). This probability is related to the probability of the previous crossing according to where we have used Eq. (42). We change variables according to χ = v 3/2 D 1/2 τ 3/2 and obtain The process is now described by a sequence of waiting times τ 1 , τ 2 , · · · and the corresponding generalized velocities v 3/2 (1), v 3/2 (2), · · ·. The displacement in the sth interval is: The advantage of this representation of the problem in terms of the pair of microscopic stochastic variables τ, v 3/2 (instead of the correlated pair τ, χ) is clear from Eq. (45): we may treat v 3/2 and τ as independent random variables whose corresponding PDFs are g(τ ) and Then P (x, t), the probability of finding the particle in (x, x+dx) at time t, is found according to Here the survival probability W (τ * ) ≡ 1 − τ 0 g(τ * )dτ * enters since the last jump event took place at t − τ * and in the time period (t − τ * , t) the particle did not cross the momentum origin. For the same reason, we have in Eq. (47) B M (.), not B(.), since the last time interval in the sequence is a meander and not an excursion. The summation in Eq. (47) is a sum of the number s of returns to the momentum origin p = 0. We note that in the analysis the particle is always moving unlike the "wait and then jump" approach used in the original CTRW model.
As usual [48,59], we consider the problem in Laplace-Fourier space where t → u and x → k. Using the convolution theorem and Eq. (45) we find This implies that reflecting the renewal property of the underlying random walk. Summing the Fourier-Laplace transform of Eq. (47), applying again the convolution theorem for Fourier and Laplace transforms and using Eq. (51), we find a Montroll-Weiss type of equation, the Fourier-Laplace transform of P (x, t): Here Ψ M (k, u) is the Fourier-Laplace transform of Eq. (52) relates statistics of velocity excursions and meanders to the Fourier-Laplace transform of the particle density. The approach is not limited to the specific problem under investigation, namely the χ ∼ τ 3/2 scaling is not a necessary condition for the validity of Eq. (52). In the general case one must revert to ψ(χ, τ ) = g(τ )p(χ|τ ) instead of the scaling form captured by B(.). Such an approach might be useful for other systems where the friction is non-linear.
At this stage, Eq. (52) still depends on since the first passage times g(τ ) from p = to the momentum origin p = 0 is dependent (see Appendix A and next Sec. for details). In fact as → 0 the number of renewals (i.e. zero crossings) tends to infinity, while in usual CTRWs, the number of renewals (or jumps) is finite for finite observation time t. In the next sections we will show how the long time results become independent of in the limit of → 0.

ASYMPTOTIC BEHAVIOR OF ψ(k, u)
The Fourier-Laplace transform of the joint distribution of jump times and lengths ψ(χ, τ ) is (53) In this section we investigate the small k and small u behaviors ofψ(k, u) which in the next sections will turn out to be important in the determination of the long time behavior of P (x, t). The small k (u) limit corresponds to large distance (time) as is well known [48].
the Laplace transform of the waiting time PDF. This waiting time PDF is investigated analytically in Appendix A. The small u behavior of g(u) differs depending on the value of D. According to Eq. (11), if D < 1 the average waiting time τ is finite and so Here τ = ∞ 0 τ g(τ )dτ is given in terms of a well known formula for the first passage time [50]. The average waiting time for a particle starting with momentum p i = to reach p = 0 for the first time is where V (p) is the effective momentum potential, Eq. (8).
Eq. (56) reflects the absorbing boundary condition at the origin and a reflecting boundary at infinity. Notice that lim →0 τ = 0 as it should and the leading order Taylor expansion yields Here Z is the normalizing partition function, Eq. (10).
Since Z diverges when D → 1 from below, we see that average waiting time diverges in that limit. For D > 1, the mean waiting time is infinite, and so Eq. (55) does not hold and instead, as we show in Appendix A,ĝ with For large τ , this yields the power-law behavior in Eq. (11), g(τ ) ∼ g * τ −(1+α) , which in fact describes the tail also for α > 1. The prefactor g * vanishes as → 0 and importantly does not depend on the full shape of the effective potential V (p), but rather only on the value of the dimensionless parameter D. The case D = 1 contains logarithmic corrections and will not be discussed here. Using Eqs. (58,59) we find non-trivial distributions of the time interval straddling time t, the backward and forward recurrence times, which were previously treated by mathematicians without the trick. This connection between renewal theory and statistics of zero crossing of the Bessel process is discussed in Appendix F.
where q(k) is the Fourier transform of the symmetric jump length distribution q(χ). Alternatively we can use Eq. (53) with u = 0, and then, upon changing variables according to χ/( √ Dτ 3/2 ) = v 3/2 and using the normalization of the scaling function B(v 3/2 ), we find Using B(k), the Fourier transform of B(v 3/2 ), Eq. (49), we can rewrite Eq. (61), This expression is the starting point for a small k expansion carried out in Appendix C which gives for ν < 2, The non-analytical character of this expansion is responsible for the anomalous diffusion of the atom's position. The first term on the right hand side is the normalization condition, the second, √ Dk ν , is consistent with the fattailed PDF of jump lengths q(χ) ∝ |χ| −(1+ν) , Eq. (11), namely an excursion length whose variance diverges since ν < 2. As expected, this behavior is in full agreement with Eq. (11) since 4/3 + β = 1 + ν. In Eq. (63), there appears the non-integer moment of the scaling function B(.), Given the formidable structure of the scaling function B(v 3/2 ), we do not describe here [67] the direct method to obtain non-integer moments like Eq. (64). Instead, we present here a method which gives v 3/2 ν indirectly.
In a future publication [67] we will discuss this and other moments of B(v 3/2 ). We can use Eq. (60) together with q(χ) ∼ q * |χ| −(1+ν) to find in Fourier spacẽ for ν < 2. In Appendix B, we investigate the area under a Bessel excursion regularized at the origin using a backward Fokker-Planck equation [6] which gives the amplitude of jump lengths Comparing with Eq. (63), we arrive at the following simple relation Thus we may use Eqs. (59,66) to find the desired νth moment of the scaling function B(.), In the limit D → ∞, we have ν = 1/3 and then We see that while the amplitudes g * and q * vanish as the convenient theoretical tool → 0, |v 3/2 | ν is independent of it in the limit, indicating the usefulness of this variable.
The second moment (v 3/2 ) 2 As we show below the second moment (v 3/2 ) 2 determines the mean square displacement of the particles. The mean vanishes since in the underlying random walk positive and negative excursions are equally likely. As for the νth moment, |v 3/2 | ν , the extraction of integer moments from the exact solution, is not straightforward. In the regime 1/5 < D < ∞, an excellent approximation is .
Similarly, for the meander we find This and Eq. (70) agrees with known results in the Brownian limit D → ∞ [8] (units and notations used in [8] are not those used by us). In Fig. 6 we show that (v 3/2 ) 2 and (v 3/2 ) 2 M , nicely match their linear approximations [67].

THE LÉVY PHASE
We now explain why Lévy statistics, and hence the generalized central limit theorem, describes the central part of the diffusion profile P (x, t) at long times for 1/5 < D < 1. The Lévy profile is cut off in the tails of the distribution, due to the correlations between jump length and time investigated here. We focus on the Lévy phase first, because it has been reported on experimentally [3].
The key idea is that in the regime D < 1, τ is finite and hence the number of jumps n scales with t/ τ when t is large [47]. At the same time the jump lengths PDF still does not have a variance (since 1/5 < D) which means that the usual Gaussian central limit theorem does not hold. Instead, due to the power-law distribution q(χ) ∝ |χ| −(1+ν) , the process belongs to the domain of attraction of a Lévy stable law. As long as x is not too large, the correlations are not important. However when x ∝ t 3/2 the simple Lévy picture breaks down, since clearly we cannot perform a jump larger than the order of t 3/2 . Thus for 1/5 < D < 1, and t −3/2 k 1, we can approximate where we have used Eqs. (55,65). Here from Eq. (65) c 0 = π/[sin(πν/2)Γ(1 + ν)]. Eq. (72) corresponds to a decoupling scheme, ψ(k, u) ĝ(u)q(k), which according to arguments in [60] is exact in the long time limit in the regime under investigation. Notice that 1/5 < D < 1 gives 2/3 < ν < 2.
Using the Montroll-Weiss type Eq.
and Ψ M (k, u) ∼ τ , which is easy to prove, we find the Fourier-Laplace representation of the solution As mentioned, both τ and q * vanish as approaches zero. Rearranging, we have where and from Eqs. (57,66), [21] K ν = lim K ν is called the anomalous diffusion coefficient. When returning to physical units, we get which has units cm ν /sec. An equivalent expression is K ν = c 0 |v 3/2 | ν D ν/2 lim →0 g * /3 τ . P (k, u) as given in Eq. (74), is in fact precisely the symmetric Lévy PDF in Laplace-Fourier space, whose (x, t) presentation (see Eq. (B17) of [62]) is for 2/3 < ν < 2. The properties of the Lévy function L ν,0 (.) are well known. The Fourier transform of this solution is exp(−K ν t|k| ν ) for 2/3 < ν < 2 which can serve as the working definition of the solution, via the inverse Fourier transform. Fig. 7 shows excellent agreement between simulations and the theory. It also illustrates the cut off on Lévy statistics which is found at distances x ∝ t 3/2 . Beyond this length scale, the density falls off rapidly. This, as noted above, is the result of the correlation between χ and τ , as there is essentially no weight associated with paths whose displacement is greater than the order of t 3/2 . This cutoff ensures the finiteness of the mean square displacement, see Fig. 8. Using the power-law tail of the Lévy PDF L ν (x) ∝ x −(1+ν) and the time scaling of the cutoff we get: for 2/3 < ν < 2 (1/5 < D < 1). If we were to rely only on the Lévy PDF, Eq. (78), we would get x 2 = ∞. Thus the Lévy PDF solution must be treated with care, realizing its limitations in the statistical description of the moments of the distribution and its tails. As we soon show, if ν > 2 we get normal diffusion x 2 ∝ t while for ν < 2/3 we have x 2 ∝ t 3 . This last behavior is due to the correlations, which restrict jumps longer than t 3/2 . So, we have three phases of motion [21]: These simple scaling arguments for the mean square displacement can be derived from a more rigorous firstprinciple approach, to which we will turn in Sec. .

The diffusion exponent
Wickenbrock, et al. [2] investigated the additional effect of a high frequency (HF) oscillating force F HF = A HF sin (ω HF t + φ 0 ) on the dynamics of the atoms, where the frequency ω HF is much larger than other frequencies in the system. According to [2] in the limit of a strong drive, the depth of the optical lattice potential is renormalized U 0 → U 0 |J 0 (2kr)| where J 0 (.) is the Bessel function of the first kind, r = A HF /(mω 2 HF ) and k is the laser field wave vector. This elegant set-up allows the control of the transport via the renormalization of the optical depth U 0 and according to Eq. (7) the control of the dimensionless parameter D. For example, in the vicinity of the zeroes of the Bessel function J 0 we clearly find an effective shallow lattice, which according to the theory corresponds to the Richardson phase. In the experiment [2] resonances in the transport are observed close to the zeros of Bessel functions, namely an enhanced spreading of the atoms. However, as pointed out in [2] many super-diffusing atoms are lost, which leads to an underestimate of the diffusion exponent. In [2] the diffusion exponent was found using Monte-Carlo simulation (see   1 there). To demonstrate the predictive power of our theory, at least for exponents, we compare between theory and the numerics [2]. As mentioned in our summary we postpone a comparison of experiments to theory until the losses become insignificant.

THE MEAN SQUARE DISPLACEMENT
Here we present the calculation of the mean square displacement of the atoms using the Montroll-Weiss equation. Our aim, as declared in the Introduction, is to unravel the quantitative connections between the transport and the statistics of excursions, for example the relation between the mean square displacement and moments of the area under the Bessel excursion and meander. Such relations are expected to be general beyond the model under investigation. A different strategy for the calculation can be based on a super-aging velocity-velocity correlation function approach [22,36]. To derive Eq. (78), we assumed in Eq. (74) that u and K ν |k| ν are of the same order of magnitude. We now consider a different small k, u limit ofP (k, u). We first expand the numerator and denominator ofP (k, u) in the small parameter k to second order. We leave u fixed and findP This second order expansion contains the second order moment of the scaling function (v 3/2 ) 2 = ∞ −∞ (v 3/2 ) 2 B(v 3/2 )dv 3/2 and similarly for (v 3/2 ) 2 M . Thus the numerator (denominator) in (81) contains a term describing the meander (excursions) contribution, respectively. From symmetry v 3/2 = 0 hence the expansion does not contain linear terms. Herê whereŴ (u) = [1 −ĝ(u)]/u is the Laplace transform of the survival probability W (τ ) and The third order derivative with respect to u is clearly related to the χ ∝ τ 3/2 scaling we have found and to the second order expansion. While the k 2 expansion in Eq. (81) works fine for small k and finite u, when u → 0 we get divergences. For example when α < 1 we haveĝ(u) ∼ 1 − Gu α + · · · and hence the third order derivative ofĝ(u) diverges as u → 0. In fact it is easy to see thatf 2 (u) will diverge when u → 0 when g(τ ) ∝ τ −(1+α) and α < 3. Thus α = 3 marks a transition from anomalous diffusion to normal which is consistent with what we found in the previous section since α = 3 gives D = 1/5 and hence ν = 2. Of course, the k 2 behavior in Eq. (83) is very different from the non-analytical |k| ν found in Eq. (63). This indicates that the order of taking the limits k → 0 and u → 0 is non-commuting [59]. The Laplace transform of the mean square displacement of the atoms is given by Hence for particles starting on the origin one can easily see The second moment is finite due to the observed fast decay of the scaling function B(v 3/2 ) for v 3/2 1 ensuring that (v 3/2 ) 2 and (v 3/2 ) 2 M are both finite.

Obukhov-Richardson diffusion
When α < 1 (D > 1), we haveĝ(u) ∼ 1 − Gu α + · · · where G = g * |Γ(−α)|. Using Eqs. (82,83),f 1 (u) ∼ c 1 u −3 andf 2 ∼ c 2 u −3 for small u with The small k, u expansion ofP (k, u), Eq. (81), iŝ which is g * and independent. The mean-square displacement in the small u limit is Converting to the time domain The scaling x 2 ∝ t 3 in this α < 1 regime is similar to Richardson's observations concerning the relative diffusion of a pair of particles in turbulence (see discussion below). We see that in this regime, both the meander and the excursion contribute to the computation of the mean square displacement. The theory agrees with the finite time simulations presented in Fig. 10, where we see Scaled mean square displacement for the Richardson-Obukhov phase. Simulations and theory nicely match without fitting. Close to the transition to the Lévy phase, i.e., D → 1, the finite time simulations slowly converge to the asymptotic limit, as might be expected. In the simulations we used t = 10 5 and averaged over 10 5 particles. that x 2 /t 3 approaches zero as D → 1 in the long time limit.
We can verify Eq. (89) in the limit of large D, where the friction force is negligible. In that case, Eq. (4) with F (p) = 0 easily gives x 2 = 2Dt 3 /3. On the other hand we have for Brownian excursions and meander (v 3/2 ) 2 = 5/6 and (v 3/2 ) 2 M = 59/30, Eqs. (70,71) [8]. Using lim D→∞ α = 1/2, we have in this limit c 1 = 15/8 and c 2 = 3/8. Plugging these numbers in Eq. (89) we get x 2 = 2Dt 3 /3, as it should. This simple demonstration, implies that it is essential to treat the last jump event properly (as a meander) and previous CTRW approaches relying on a unique jump length distribution, lead to wrong conclusions (e.g. replacing (v 3/2 ) 2 M with (v 3/2 ) 2 is wrong). Furthermore, in this limit the contribution of the meander is numerically larger then the contribution of the excursions, even though the number of excursions is large. Notice however that for the calculation of the Lévy density P (x, t), Eq. (78), the statistics of the meander did not enter. Hence depending on the observable of interest, and the value of D the meander may be either a relevant part of the theory or not.

Super-diffusion
For 1 < α < 2 (1/3 < D < 1) we use the expansion Here the first moment of the waiting time is finite while the second moment diverges. Crucially both τ and G = |g * Γ(−α)| vanish as → 0. For this parameter regime, we have found Lévy behavior for the central part of P (x, t). Using Eq. (85) we find (91) Using Eqs. (57,59) The inversion of Eq. (91) to the time domain and inserting D = (2α − 1) −1 gives (93) The same result is valid for 2 < α < 3. This behavior depends on the normalizing partition function Z, namely the shape of the potential V (p) in the vicinity of the momentum origin becomes important, unlike the Richardson-Obukhov phase, where only the large p behavior of V (p) is important in the long time limit (i.e., the case α < 1, Eq. (89)). Notice that |c 1 | in Eq. (93) tends to zero when α → 3 from below. This means that the meander becomes irrelevant when we approach the normal diffusion phase x 2 ∼ t. Fig. 11 compares simulation and theory and demonstrates that x 2 /t 4−α diverges as the transition to the Gaussian phase D < 1/5 is approached.

Breakdown of scaling assumption-Normal Diffusion
Our starting point was the scaling hypothesis χ 2 ∼ τ 3 , Eq. (14) and Fig. 3. Indeed we have shown that for large χ and τ the scaling function B(.) describes the conditional probability density p(χ|τ ) for all values of D. However this does not imply that the scaling solution B(.) is always relevant for the calculation of the particle density P (x, t). Roughly speaking, so far we have assumed that large jumps χ and long waiting times τ dominate the underlying process. No one guarantees, however, that this large χ limit (long jumps) is important while small χ (small jumps) can be neglected. Indeed when we switch over to the normal diffusion phase D < 1/5, the small scale aspects of the process become important (meaning the regularization of the potential V (p) when p → 0 becomes crucial). This is similar to the Lévy versus Gaussian central limit theorem arguments, where the former is controlled by the tails of the distribution while the latter by the variance. Let us see this breakdown of scaling in more detail.
When α > 3 we have the expansionĝ(u) ∼ 1 − u τ + u 2 τ 2 /2−u 3 τ 3 /6+· · · where the first three integer moments of the waiting time PDF are finite (see Appendix A). Then the functionf 1 (u) is negligible when u → 0 whilef 2 (u) ∼ τ 3 /u τ . We find Thus α = 3 (or D = 1/5) marks the transition between the anomalous super-diffusive phase and the normal diffusion phase. Furthermore, in this case the meander is of no importance. However, Eq. (94) is correct only if our assumptions about scaling are valid. By definition and if the scaling hypothesis holds which gives Inserting Eq. (97) in (94) we get the expected result, discussed in the introduction This simple result is the correct one, even though the scaling assumption is not valid in this regime. Sometimes we reach truthful conclusions, even though the assumptions on the way are invalid. To see the breakdown of scaling we plot in Fig. 12 Y dev = lim →0 χ 2 /[D (v 3/2 ) 2 τ 3 ]−1 versus 1/D which should be zero in the normal phase if the scaling hypothesis holds (where χ 2 and τ 3 are finite). The calculations of χ 2 and τ 3 are given in Appendix A and B. We see that at the transition point, D = 1/5, Y dev = 0, but otherwise Y dev = 0. Hence the scaling hypothesis does not work in the normal phase D < 1/5. This means we need another approach for normal diffusion, which luckily is easy to handle. For D > 1/5 long jumps dominate, since χ 2 diverges, and then our scaling theory works fine.

THE RICHARDSON-OBUKHOV PHASE
In Sec. we discussed the Lévy phase which is found for 1/5 < D < 1. When the average jump duration, τ , diverges, i.e. for D > 1, the dynamics of P (x, t) enters a new phase. Since the index of the Lévy PDF ν approaches 2/3 as D approaches 1, x scales like t 3/2 in the limit. Due to the correlations between χ and τ , x cannot grow faster than this, so in this regime, P (x, t) ∼ t −3/2 h(x/t 3/2 ). An example for this behavior is presented in Fig. 13. This scaling is that of free diffusion, namely momentum scales like p ∼ t 1/2 and hence the time integral over the momentum scales like x ∼ t 3/2 . Indeed in the absence of the logarithmic potential, namely in the limit D 1 the Langevin equations (4) give  99) is valid when the optical potential depth is small since D → ∞ when U 0 → 0. This limit should be taken with care, as the observation time must be made large before considering the limit of weak potential. In the opposite scenario, U 0 → 0 before t → ∞, we expect ballistic motion |x| ∼ t, since then the optical lattice has not had time to make itself felt [3]. Physically, the atoms in this phase are undergoing a random walk in momentum space, due to random emission events, which in turn give the x ∼ t 3/2 scaling. For shallow lattices, the Sisyphus cooling mechanism breaks down, in the sense that transitions from maximum to minimum of the potential field created by the laser fields, are not preferred over the inverse transitions. Thus the deterministic dissipation is not effective, and we are left with Brownian scaling in momentum space, p ∼ t 1/2 .

THE NORMAL PHASE
When the variance of the jump length is finite, namely ν > 2 (D < 1/5), we get normal diffusion. Here the variance of jump lengths is finite and hence the scale free dynamics breaks down. The breakdown of scaling means that instead of using the scaling function B(.), e.g. in Eq. (42) for ψ(χ, τ ) we must use the joint PDF ψ(χ, τ ) = g(τ )p(χ|τ ) Eq. (13) and in principle not limit ourselves to large τ . However, luckily there is no need for a new calculation.
We focus on the central part of of the density P (x, t) where central limit theorem arguments hold. In this normal case the spatio-temporal distribution of jump times and jump lengths effectively decouples, similar to the Lévy phase. Since the variance of jumps size and the averaged time for jumps are finite, many small jumps contribute to the total displacement, and hence in the long time limit, we expect Gaussian behavior with no correlations between jump lengths and waiting times, i.e., the decoupling approximation is expected to work. More precisely, the average waiting time is finite sô g(u) ∼ 1 − τ u + · · · and the variance of jump lengths is also finite so the Fourier transform of q(χ) has the following small k expansionq(k) = 1 − χ 2 k 2 /2 + · · · where χ 2 = ∞ −∞ χ 2 q(χ)dχ is the variance of the jump lengths. This variance is investigated in Appendix B using a backward Fokker-Planck equation. In the small k, u limitψ(k, u) ∼ 1 − u τ + k 2 χ 2 /2 + · · · and the Montroll-Weiss equation (52) iŝ This is the expected Gaussian behavior for the position probability density and is the diffusion constant, namely Eq. (102) relates the statistics of the excursions, i.e. the variance of the area under the excursion χ and its average duration τ to the diffusion constant (here χ = 0). As noted in the introduction the equation has the structure of the famous Einstein relation, relating the variance of jump size and the time between jumps to the diffusion. While τ , Eq. (57), and χ 2 Eq. (155), approach zero when → 0, their ratio remains finite and gives (104) This equation was derived previously using different approaches [20,36]. The diffusion constant K 2 diverges as D → 1/5 from below indicating the transition to the super-diffusive phase. A sharp increase in the diffusion constant K 2 as the intensity of the laser reaches a critical value was demonstrated experimentally [20]. Eq. (104) can be derived using the Green-Kubo formalism [6], so in the normal phase, the analysis of the statistics of excursions is an alternative to usual methods. It seems that in the Lévy phase, the analysis of statistics of excursions is vital. Specifically the usual Green-Kubo formalism breaks down since K 2 is infinite, and the calculation of the anomalous diffusion coefficient K ν cannot be based on a computation of a stationary velocity correlation function.

DISCUSSION
Starting with the Langevin description of the semiclassical motion of atoms undergoing a Sisyphus cooling process, we mapped the problem onto a random walk scheme using excursions as a tool. Thus our work combines ideas from the theory of stochastic processes with cold atoms physics. We now summarize the key ingredients of our results and its predictions from the point of view of these two communities.

Fractional Diffusion Equation and CTRW theory revisited
The first ingredient of the theory was the calculation of the coupled joint distribution ψ(χ, τ ). Rather generally, coupled waiting time -jump length distributions are the microscopical ingredient for coupled continuous time random walks [48,[58][59][60][61][62][63]. These distributions, however are difficult to obtain from first principle calculations and usually treated in a simplified manner. For example, postulating ψ(χ, τ ) = g(τ )δ(|χ|−τ )/2 is common in stochastic theories. Hence we provided an important pillar for the foundation of this widely applied approach.
Our work gives the relation between velocity excursions, and random walk theory and hence diffusion phenomenon. Since zero crossing in velocity space is obviously a very general feature of physical paths of random processes, we expect the areal distributions of excursions and meanders will play an important role in other systems, at least as a tool for the calculation of transport and diffusion (e.g., in principle our approach can treat also the case of a constant applied force, where a mean net flow is induced). The celebrated Montroll-Weiss equation needed two important modifications. First, since the underlying process is continuous we regularized the process, such that the excursions and meander start at → 0. Secondly, the statistics of the last jump event, i.e., the meander, must be treated with care. In contrast, usually it is assumed that only ψ(χ, τ ) is needed for a microscopical description of a continuous time random walk model. This and the correlations between χ and τ make the problem challenging and usual approaches to diffusion fail. With our approach, the behavior of the packet is mapped onto a problem of the calculation of areas under Bessel excursions and meanders. For example we derived the relation between the mean square displacement of the atoms and the areas under both the Bessel excursion and meander Eqs. (85, 89, 93). A decoupled scheme which neglects the correlations ψ(χ, τ ) g(τ )q(χ) gives a diverging result for ν < 2, which is unphysical.
In the regime 1/5 < D < 1, which we have called the Lévy regime, the correlations between jump lengths and waiting times point to the limitations of the fractional diffusion equation, a popular framework based on frac-tional calculus [48]. The fractional diffusion equation was previously investigated in the context of random walk theory, and it describes a Lévy flight processes [48]. Here we have provided a microscopic justification for it. The fractional diffusion equation [48,71,72], was the phenomenological starting point for the description of the experiments in the work of Sagi et al. [3].
Our work (see also [21]) provides the exponent ν, Eqs. (7,30), in terms of the recoil energy and lattice depth, β = 1 and the anomalous diffusion coefficient, K ν from Eq. (76). Indeed, the experiment [3] found the value β = 1, so the time derivative on the left hand side is a first order time derivative. The fractional space derivative ∇ ν is a Weyl-Reitz fractional derivative [48]. To see the connection between our results and the fractional equation we re-express Eq. (74) as which is the Fourier-Laplace transform of Eq. (105).
Here we recall [48] that the Fourier space representation of, ∇ ν , is −|k| ν and that u P (k, u) − 1 is the representation of the time derivative in Laplace-Fourier space, of a δ function initial condition centered on the origin. We see that P (x, t) for 2/3 < ν < 2 satisfies the fractional diffusion equation. In other words, for initial conditions starting on the origin the solution of Eq. (105) is the Lévy PDF Eq. (78). However the use of the fractional diffusion equation must be performed with care. It predicts a seemingly unphysical behavior: the mean square displacement is infinite. Indeed as mentioned in the introduction, the mean square displacement was shown experimentally to exhibit superdiffusion by Katori et al. [1] and in simulations in [2]. In fact the mentioned coupling implies that the fractional equation is valid only in a scaling region, |x| < t 3/2 and thus the Lévy distributions obtained analytically here and used phenomenologically in the Weizmann experiment [3] describing the center part of the spreading distribution of the particles, but the tails exhibit a cutoff for x > t 3/2 . In other words the Lévy distribution describing the center part of the packet for 1/5 < D < 1 does not contain information on the correlations and to experimentally investigate correlations in this regime one must probe the tails of the packet.
The renewal approach is the basis of the coupled CTRW theory developed here. Renewal theory turned out to be predictive in the sense that while the set of zero crossings is nontrivial in the continuum limit, we could avoid this mathematical obstacle by introducing modified paths with an ± starting point after each zero hitting. Physical observables, like ν, P (x, t), K ν and x 2 do not depend on in the limit, as expected. Technically, this is due to a cancellation of , for example both τ , Eq. (57), and q * , Eq. (66), depend linearly on hence their ratio is independent and the anomalous diffusion coefficient K ν , Eq. (76), becomes insensitive to . The same is true for the normal diffusion constant K 2 Eq. (102) and other physical observables Eqs. (69,89,92,93). A closer reading of Appendix A and B will reveal that this cancellation is not exactly trivial. In our work we use the regularized Bessel process for the calculation of the PDFs of the first passage time, g(τ ), and jump lengths, q(χ). Hence for the calculations the details of the potential V (p) for small p are important (but not for B(v 3/2 )!). If we use instead the PDFs g B (τ ) and q B (χ) (see Appendix A and B) of the non-regularized process, i.e. a logarithmic potential for any p, some (but not all) of the cancelation will not take place (see Eq. (67) where both q * ∝ and g * ∝ for the regularized process which allows for the cancelation). For example both for the regularized and the non-regularized processes g(τ ) ∼ g * τ −(1+α) (so the exponent α is identical in both cases), however g * ∝ for the regularized process, while g * ∝ 2α for the nonregularized case Eq. (113). The cancelation of is further discussed in Appendix F for three types of random durations. That Appendix shows that while the first passage time PDF g(τ ) depends on (see Appendix A), the time interval straddling time t and the backward and forward recurrence times are insensitive to this parameter. Thus one upshot of Appendix F is a further demonstration that this intuitive renewal trick actually works. A rigorous mathematical treatment of the problem will lead to a stronger foundation of our results.
Regularization of the process is crucially important in the Lévy and Gaussian phase D < 1 where observables depend on the details of the potential V (p). This is related to work of Martin et al. [41] on the classification of boundaries for the non-regularized Bessel process. They showed that, using our notation, for D < 1, p = 0 is an exit boundary while for D > 1 it is a regular boundary (D < 0 corresponds to an entrance boundary condition, which is not relevant to our work). For an exit boundary starting on p 0 it is impossible to reach a finite momentum state p so clearly we cannot afford such a boundary in our physical problem. For that reason, we need to consider the regularized first passage time problem at least when D < 1, and then the boundary is regular. A regular boundary on p = 0 means the the diffusion process can enter and leave from the boundary. Therefore, for D > 1 our final results do not depend on the shape of V (p), besides its asymptotic limit, namely one may replace the regularized Bessel process with the non-regularized one. To see this, notice that x 2 for D > 1 depends on α < 1 but not on small scale properties of V (p) and since α is the first passage exponent for both the regularized and non-regularized processes, it does not really matter which process is the starting point of the calculation. Further evidence for such a behavior is in Appendix F, where statistics of durations for D > 1 are shown to depend on α < 1 but not on the shape of V (p). Thus the D = 1 marks transition from Richardson to Lévy behavior, the diverging of the partition function Z namely the equilibrium velocity distribution turns non-normalizable, and according to Martin et al. the boundary of the Bessel process switches from regular to exit.
Additional theoretical work is needed on the leakage of particles, and more generally evaporation (see more details below). Assuming that energetic particles are those which get evaporated, and that this takes place through the boundary of the system, one would be interested in the first passage time properties of particles, from the center of the system to one of its boundaries. Previous theoretical work on first passage times for Lévy walks and flights might give a useful first insight, [73][74][75][76][77] however it is left for future work to compare these simplified models with the microscopic semiclassical picture of the underlying dynamics. Evaporation is an important ingredient of cooling (beyond the Sisyphus cooling) since it gets rid of the very energetic particles; hence this line of research has practical applications. Theoretically, understanding the boundary conditions for fractional diffusion equations, needed for the calculation of statistics of first passage times, is still an open challenge. Yet another challenge is to characterize the joint PDF W (x, p, t) and in particular the correlations between position and momentum, which are expected to be non-trivial. Again, some elegant ideas and tools were recently developed [78] within the Lévy walk framework, however more work is needed for direct comparison with the physical picture under consideration.

Cold atom experiments
Experiments on ultra-cold diffusing particles in optical lattices have been performed on a single 24 Mg + ion [1] and an ensemble of Rb atoms [2,3,20]. Experiments on the spreading of the density of atoms yield ensemble averages like the mean square displacement and the density P (x, t). Analysis of individual trajectories yields, in principle, deep insight on paths, for example information on the first passage time in velocity space, i.e. g(τ ), zero crossing times, statistics of excursions and meanders etc. For a comprehensive understanding of the dynamics of the particles, both types of experiments are needed. The theory presented herein provides statistical information on the stochastic trajectories, e.g., zero crossing events, which in turn are related to the spreading of the packet of particles. Thus from single trajectories, one may estimate g(τ ), B(v 3/2 ), and ψ(χ, τ ). From g(τ ) one may then obtain the exponent α which can be used to predict the qualitative features of the spreading of the ensemble of particles. For example with an estimate of α we can determine the phase of motion, be it Richard-son, Lévy or Gauss. Of course, a more quantitative investigation is now possible, since we have analytically related α, the mean square displacement, and K ν with microscopical parameters like the optical potential depth U 0 . In principle, statistics of Bessel excursions, so far of great interest mainly in the mathematical literature, could be detected in single particle experiments. These single ion experiments are also ideal in the sense that collisions/interparticle interactions, do not play a role and they also provide insights on ergodicity. Three specific unsolved issues are: (i) Can the Richardson phase be experimentally demonstrated? So far this phase was obtained in our theory, and Monte Carlo simulations [2], while experiments exhibit ballistic diffusion as the upper limit [1][2][3]. This might be related to the subtle limit of taking time to infinity before the depth of optical lattice approaches zero, since clearly in the absence of an optical lattice the fastest particles are ballistic.
(ii) While our theory predicts correctly a Lévy phase in agreement with experiments, in [3] two exponents were used to fit the data. In contrast the theory we developed suggests a single Lévy exponent ν. This could be due to the leakage of particles and also to the t 3/2 cutoff of the Lévy density, which are due to the correlations between χ and τ investigated in this manuscript. Avoiding the leakage of particles is an experimental challenge, and once this is accomplished, a more informed comparison of our theoretical prediction on the Lévy phase with a single characteristic exponent ν could be made. Success on this front would provide an elegant demonstration of Lévy's central limit theorem, with the characteristic exponent controlled by the depth of the optical potential.
(iii) As pointed out in [2], losses of atoms lead to an underestimate of the diffusion exponent, as many super-diffusing atoms are lost. Characterisation of these losses both theoretically and experimentally could advance the field, since this yields insight on the underlying processes and leads as towards better control of the particles. Specifically, we do not know what is the number of particles kept in the system, versus time, for varying strength of the optical potential.
(iv) In the Lévy phase, measurements of non-integer moments |x| q with q < ν < 2 will according to theory exhibit Lévy scaling x 2 ∼ t q/ν . In contrast, the main focus of experiments so far was the second moment q = 2, which is difficult to determine statistically since as we have shown here it depends on the tails of P (x, t) and very fast particles. Thus ideally a wide spectrum of moments |x| q should be recorded in experiment, the low order moments q < ν giving information on the center part of the density while higher order moments on the correlations and tails.
In the case that experiment and theory will not reconcile, we will have strong indication that the current semi-classical theory in not sufficient and then we will be forced to investigate at least four other aspects of the problem: (a) Effects of collisions on anomalous diffusion of the atoms.
(b) Effect of higher dimensions.
(c) Quantum effects beyond the semi-classical approach used here. In particular it would be very interesting to simulate this system with full quantum Monte Carlo simulations [6], to compare the semiclassical theory with quantum dynamics. We note that the Richardson phase, which as mentioned was not observed in experiments, is actually a heating phase and quantum simulations become difficult because the numerical lattice introduces a cutoff on velocities which induces artificial ballistic motion.
(d) Other cutoff effects that modify the anomalous diffusion. For example at high enough velocities, Doppler cooling is expected to diminish the fast particles. In [79] we simulated the effect of Doppler friction, and showed that the anomalous character of the diffusion is kept unchanged, at least for a certain reasonable set of parameters. However, a general rule on the influence of Doppler cooling is not yet established, and an experimentalist with a specific set of parameters in mind might wish to test numerically the magnitude of this effect on the anomalous spreading.
An interesting approach was recently suggested by Dechant and Lutz [23]. They investigated the multifractal nature of the moments |x| |q| of the process and consider initial conditions different than ours. They assume that in a first stage, of duration t c , the particles are cooled in a confining field (which inhibits the spreading). Then the momentum of particles relaxes to a state described by the infinite covariant density [26] which depends on t c . The particles are then released and their spreading is recorded for a duration t. We considered the case t t c while Ref. [23] considers the opposite case. At-least the experiment in [3] is conducted under the conditions investigated here, namely that the spreading time is much longer than the preparation time. As shown by Hirschberg et al., starting with power law distributions will dramatically influence the spreading, both in momentum space [25] and in space. Indeed large t c implies power law initial conditions, Eq. (9), with a t c dependent cutoff [26]. In that sense diffusivity is sensitive to the initial preparation of the system, and an experimental verification of these effects would indicate the fundamental difference between transport in these systems and normal transport which does not depend on initial conditions. This is clearly related to the strong sensitivity we have found of the mean squared displacement on a single jump event, described by the Bessel meander.
Acknowledgement This work was supported by the Israel Science Foundation. We thank Yoav Sagi and Nir Davidson for discussions on the experiment [3] and Andreas Dechant and Eric Lutz on collaborations on related problems.

The waiting time PDF
The waiting time is the time it takes the particle starting with momentum p i > 0 to reach p f < p i for the first time. Here we investigate its PDF, g(τ ).

First passage time for the Bessel process
We first briefly investigate the first passage time problem for the Bessel process following Ref. [35]. The Bessel process corresponds to the case F (p) = −1/p so the the force diverges at the origin. In the next subsection we will consider the regularized force F (p) = −p/(1 + p 2 ). According to [50], the survival probability S(τ ), namely the probability that a particle initially at p i does not cross the boundary p f < p i in the time interval whose length is τ , satisfies Here S = 1 for τ = 0 since the particle's escape time cannot be zero. Further S → 0 when p i → p f , since in that case the particle starts out at the boundary, and S = 1 if one starts at p i = ∞. The random time it takes a particle starting on p i to reach p f for the first time is τ and its PDF is g B (τ ) = −∂ τ S(τ ). Here the subscript B denotes the Bessel process. Since S(τ )| τ =0 = 1 we have in Laplace τ → u space the following simple relation Using Eq. (108) and the Laplace transform of Eq. (107) we find The solution of Eq. (109) is [80] g B (u) = Here K α (.) is the modified Bessel function of the second kind and as before α = 1/2 + 1/(2D). If α > 1 the small u expansion of Eq. (110) yieldŝ g B (u) ∼ 1 − u τ B + · · · where τ B is the average first passage time. Expanding Eq. (110) we find Notice that τ B diverges when α → 1 corresponding to D → 1. Not surprisingly, the average time in Eq. (111) is generally different from the expression for the average waiting time for the regularized process Eqs. (56,57). Expanding Eq. (110) for small u for the case 0 < α < 1 we get (112) Inverting to the time domain we find, using a Tauberian theorem, the long time behavior of the PDF, One can show that this fat-tail behavior is valid also for the regime α > 1. The procedure involves expansion of Eq. (110) beyond the first two leading terms; for example for 1 < α < 2 one findsĝ B (u) = 1 − u τ B + Cu α .... and the third term gives the tail of PDF. Note that the PDF g B (τ ) yields the same exponent as in Eq. (11). However to calculate the amplitude g * we must consider the regularized process.

First passage time for the regularized process
For the optical lattice problem we need to treat the regularized force F (p) = −p/(1 + p 2 ). Our aim is to find the asymptotic behavior in Eq. (11) while the average waiting time is given already in Eqs. (56,57). From these, we know that for α < 1 the average first passage time from p i to p f is infinite (the derivation can be easily generalized for arbitrary initial and final states). Furthermore, for long time the large momentum behavior of F (p) plays the crucial role and hence we expect for small uĝ where we must determine G which then gives the amplitude g * introduced in Eq. (11). The equation forĝ(u) is the same as Eq. (109) but with the regularized force: We need to solve forĝ(u) for small u, so to leading order we drop the last term and find where V (p) = ln(1+p 2 ) 1/2 . Sinceĝ(u) = 1 when p i = p f , C 2 (u) = 1. This approximation breaks down however when p i is too large, since then,ĝ(u) ∝ p 1+1/D i , and the last term in Eq. (115) is comparable in size to the first two when p i ∼ u −1/2 . In the large p i regime, however, p i 1, so we can approximate F (p i ) by its Bessel form, F (p i ) ≈ −1/p i , and the solution is, as above, where here we cannot use the p f → p i limit to normalizeĝ(u), since p f is not necessarily large. Nevertheless, g(u = 0) = 1 and K α (z) ∼ (Γ(α)/2)(z/2) −α for z → 0, implying that to leading order in u. For 1 p i u −1/2 , our two approximations must agree, and so we find using the second order expansion of K α (z) and recognizing that (1/D + 1) = 2α, Thus, using Eq. (117), α pi p f e V (p )/D dp u α .
(122) Using the Tauberian theorem, which implies u α → τ −(1+α) /Γ(−α), we find the large τ behavior It can be shown, with additional work, that this result is also valid for α > 1. To conclude, we find for the final state p f = 0 and the initial state p i = 1, This of course vanishes as → 0, in this case linearly, as opposed to the p i → and p i = 0 limit of Eq. (113), where the dependence is of higher order. As shown in the manuscript, for the purpose of calculation of the asymptotic behavior of P (x, t), all we need is ∂ g * when → 0 which is a finite constant which depends on D only. The full shape of V (p) is unimportant for the calculation of g * indicating a degree of universality. In contrast τ depends on Z and so is non-universal. Since g * is a measure of the long time behavior, the detailed shape of the potential is not important, provided it is regularized.

Moments of g(τ )
In the text we show that τ and τ 3 are relevant when D < 1/5 so that these averages do not diverge. The moments τ and τ 3 are found usinĝ is the Laplace transform of the waiting time PDF g(τ ). From Solving Eq. (127) so τ = 0 when p = 0 as expected for the first passage from p to the origin. For the second and third moments and Thus, in the p = → 0 limit, The jump length PDF q(χ) The excursion length χ is the distance the particle travels from its start with momentum p i > 0 until it reaches the momentum origin p f = 0 for the first time. Here we investigate its PDF q(χ) which in the limit of p i → 0 and for the regularized force field gives the jump length PDF. Previously this PDF was investigated for the −1/p force field in [6]. Note that in what follows p i > 0 so χ > 0. As elsewhere in this paper, from symmetry positive and negative χ are equally probable so we may restrict our attention to χ > 0.

PDF of jump lengths for the Bessel process
As mentioned, for the Bessel process F (p) = −1/p. The PDF q B (χ) which depends of course on the start and end points p i and p f satisfies the following backward equation [50]: where the subscript B again stands for Bessel. Since χ > 0, we define the Laplace transform When p f → p i , we have q B (χ) → δ(χ) and q B (χ) χ=0 = 0 if p i = p f since it takes time for the particle to reach the boundary and hence we cannot get an excursion whose size is zero. Of course, in the opposite limit of large χ, lim χ→∞ q B (χ) = 0. In Laplace space Eq. (134) yields It is easy to verify that the appropriate solution iŝ with ν = (1 + D)/3D and as in the previous Appendix K ν (·) is the modified Bessel function of the second kind.
Our goal is to find the large χ behavior of q B (χ), so we expand Eq. (137) in the small s limit, finding, for 0 < ν < 1, The first term on the left hand side is simply the normalization, and the s ν term indicates that the average excursion length diverges since ν < 1. Passing from s to χ, we get for large χ With a similar method, one can show that Eq. (139) holds also for 1 < ν < 2. There the expansion Eq. (138) contains three terms; the additional term yields the average of χ which is now finite. In [6] a more complicated method was used to investigate the same problem and a similar result was found although with a typographical error. They report (3ν + 1) ν in the prefactor whereas we find (3ν − 1) ν , a difference with some importance since the pre-factor found here goes to zero when D → ∞, which is necessary to obtain reasonable physical results.

PDF of χ for regularized process
To obtain the large χ behavior of q(χ) for the regularized force F (p) = −p/(1 + p 2 ), we follow the same steps carried out in Sec. . The equation to solve is For ν < 1, we switch to Laplace space χ → s, i.e. D(∂ pi ) 2q + F (p i )∂ piq − p i sq = 0, and drop the last term for small s, yieldingq after applying the boundary conditionq| p f →pi = 1. Again, for p i 1, we can approximate F (p i ) by its large p i approximation, yieldingq In the small s limit, q(s → 0) = 1 implies to leading order. These two approximations must agree for 1 This implies that This calculation was done assuming p i , p f positive. Allowing also for negative momenta, and so negative χ, we get an additional prefactor of 1/2: One can show that this result is valid in the whole domain of interest 1/3 < ν < 2, i.e. the domain where the variance of χ is infinite.
The variance of χ We here obtain the finite χ 2 for the case ν > 2 for the regularized force F (p) = −p/(1 + p 2 ). Using the Laplace χ → s transform of Eq. (140) we find the following backward equation forq(s): Here we used q(χ)| χ=0 = 0 since the particle starting with p i > 0 cannot reach zero momentum without traveling some finite distance. The Laplace transformq(s) is expanded in s: Here χ is the average jump size, for positive excursions i.e. those that start with p i > 0 and χ 2 is the second moment. Note that in the original model we have both positive and negative excursions and so we have p i = + or p i = − with probability 1/2 hence for that case −∞ < χ < ∞ and from symmetry χ = 0. The variance of the original process is χ 2 due to symmetry. As we restrict ourselves to p i > 0, χ > 0 and χ in Eq. (148) is finite. Inserting Eq. (148) in Eq. (147) we find The s 1 terms give while the s 2 terms are D 2 This means that the first and then the second moment can be found by repeated integration. The boundary conditions are that both the first and second moments are zero if p i = 0 while they diverge if p i = ∞. Using a reflecting boundary at p i = ∞ we find (see e.g. [50], chapter XX) where the lower bound is the momentum origin and the potential is given in Eq. (8). Integrating Eq. (151) we find As explained in the text we consider only the limit when p i = is small so we have since V (0) = 0. Using Eq. (152) and integrating by parts we get Note that from symmetry V (p) = V (−p), so we have Since g(τ ) ∼ g * τ −3/2−γ for large τ , with γ = 1/(2D), we assume that there exists some t 0 such that g(τ ) = g * τ −(3/2+γ) for τ > t 0 . Then using Eq. (59) and ν = 1/3 + 2γ/3, The first term on the right hand side yields the normalization condition ψ(k, u)| k=u=0 = 1, the second term is zero for k = 0 and due to symmetry of the jumps B(v 3/2 ) = B(−v 3/2 ), this term when expanded in k gives a k 2 term. As we will show, the last term gives a |k| ν term, so as long as ν < 2 we can neglect the second term. Then changing variables tok = k √ Dτ 3 we find Here we have taken the small k limit in such a way that k(t 0 ) 3/2 which appears in the lower limit of an integral is negligible. From Eq. (157) we must investigate the integral From symmetry B(v 3/2 ) = B(−v 3/2 ) and hence the Fourier transformB(k) is a real function so and k c will eventually be taken to zero. Using the definition of the Fourier transform Let iv 3/2 = −x and then using integration by parts, Inserting the last expression in Eq. (160), the diverging (k c ) −ν terms cancel each other, and we find Using x = −iv 3/2 we get where v 3/2 ν = ∞ −∞ |v 3/2 | ν B(v 3/2 )dv 3/2 . Inserting Eq. (163) in (157), using (158) we get Eq. (63).

On the simulations
The simulation In this appendix, we briefly discuss the simulations presented in the paper to test our predictions. We have two classes of simulations: one generates Bessel excursions, and the second solves the Langevin equations. The first is based on a discrete random walk treatment of the excursion, wherein the particle takes a biased walk to the right or left in momentum space every time step. The de-gree of bias varies with p in accord with the force F (p). Our dimensionless control parameter D is the ratio of the coefficient of the diffusion term to the coefficient of the 1/p large p behavior of the force. Since for our simple random walk the diffusion constant is 1/2, the parameter D enters via the strength of the bias. The probability to move right or left is then given by The continuum limit is approached, as usual, when ∆ p → 0. Given that in our units F (p) varies on the scale of unity, it is sufficient to take ∆ p 1. In practice, we took ∆ p = 0.1. We also monitor the position, adding p to the position every time the particle takes a step starting at momentum p.
To simulate paths of the Langevin dynamics we used both a straightforward Euler-Maruyama integration of the Langevin equation as well as one based on the random walk picture above. Both methods give equivalent results.

The Bessel Excursion
An efficient way to generate the set of Bessel excursions is to generate all discrete Brownian excursions with equal probability and then to weight each excursion by its appropriate weight. The task of generating the ensemble of Brownian excursions is at first sight a nontrivial problem, since the constraint that the random walk not cross the origin is nonlocal in the individual left/right steps constituting the walk. However, it is easily accomplished using Callan's proof [81] of the famed Chung-Feller theorem [82]. Callan gives a explicit mapping of any N left, N right step random walk to a unique walk that does not cross the origin. This mapping maps exactly N + 1 of the 2N N N left, N right walks, to each of the nonzero crossing walks (of which there are 2N N /(N + 1)). Since every Brownian excursion of N + 2 steps consists of an initial right step, a non-zero-crossing N left, N right step walk and a final left step back to the origin, Callan's mapping allows us to generate an N + 2 step Brownian excursion by generating a random N left, N right step walk, apply the map, and pre-and post-pending the appropriate steps. In this way every excursion is generated with exactly the same probability.
Callan's mapping is as follows. For a given N left, N right step walk, one accepts it as is if it does not cross the origin. Otherwise, one finds the leftmost point reached by the walk (i.e., the minimal x of the path), which we denote x n , where n is the number of steps to reach this point. If this point is visited more than once, then one takes the first visit. One then constructs a new walk based on the original walk. We divide the original walk  Table I. The simulation method is outlined in Appendix and we averaged over 2×10 5 samples with τ = 10 5 .
Theoretical curve was plotted with Maple.
into two segments, the first n−1 steps and the remaining part. The first part of the new walk consists of the second segment, shifted by −x n such that the new walk starts on the origin. In the original walk the step before reaching x n is obviously a left step. One appends this left step at the end of the walk under construction, and then attaches the first segment (of the original walk). This new walk clearly does not cross the origin, while preserving the number of left and right steps, hence it starts and ends on the origin. To account for the bias, i.e., to switch from Brownian to Bessel excursions, it is sufficient to weight each excursion by the weight factor, F : Step i is a right step 1 −

F (pi)∆p 2D
Step i is a left step (165) This works well as long as D is not too small. For very small D, the convergence is quite slow since small areas have anomalously large weight unless N is very large. As long as one is interested in the large N behavior, better convergence is achieved by taking F (p) = θ(p − 1)/p.

Areal distribution of the Bessel Meander
Our goal is to find the conditional PDF p M (χ * |τ * ) where χ * = τ * 0 p(t )dt , is the area under the Brownian meander. The starting point for the calculation of p M (χ * |τ * ) is a modification of Eq. (20). As for the area under the Bessel excursion, we consider only the positive meander where 0 < χ * < ∞, and later use the symmetry of the process to find the areas under both positive and negative meanders. Since contrary to the condition on the excursion, the meander is not bound to return to the origin, we now keep the end point free and integrate the propagatorĜ τ * (s, p|p i ) over all possible values of p. Hencep M (s|τ * ) = lim →0 ∞ 0Ĝ τ * (s, p| )dp ∞ 0Ĝ τ * (s = 0, p| )dp ; (166) We expand the propagatorĜ τ * (s, p|p i ) in a complete orthonormal basis, using the same approach as in the main text, while accounting for the boundary conditions of the meander. We thus rewrite Eq. (20) for the case of the meander and integrate over p ∞ 0Ĝ τ * (s, p| )dp = 1/2D s 1/3 where the functions f k (.), as before, are the solutions of the time independent equation (27) and E k = s 2/3 λ k . Changing variables top = s 1/3 p, and using the small p behavior Eq. (28), which gives for the numerator in Eq. (166). We now define a new k-dependent coefficient The a k 's, similar to thed k 's, are evaluated from a numerically exact solution of Eq. (27). Rewriting in terms of the new coefficients ∞ 0Ĝ τ * (s, p| )dp = 2α s ν kd k a k e −Ds 2/3 λ k τ * .

27(χ
This is the main result of this Appendix. In Fig. 14) we plot this areal distribution, comparing it with a histogram obtained from finite time simulations.
The time interval straddling t The velocity process v(t ) restricted to the zero-free interval containing a fixed observation time t is called the excursion process straddling t, and the portion of it up to t is the meandering process ending at t. For Brownian and Bessel processes these random paths were the subject of intense mathematical investigation, e.g. [17,31,83] and references therein. The investigation of the duration of the excursion straddling t, the time interval for the meander ending at t, and statistical properties of the path e.g., its maximal height, have attracted mathematical attention since these reveal deep and beautiful properties of Brownian motion in dimension d (as is well known the Bessel process is constructed from the radius of a d dimensional Brownian motion). One aspect of the problem is the quantification of the properties of the set of points on the time axis on which the zero crossings take place. Clearly the number of zero crossings for a Brownian path starting at the origin within the time interval (0, t) is infinite, due to the continuous nature of the path (see more details below). One might naively expect that the time between points of zero crossing approaches zero, since a finite measurement time divided by infinity is zero. That, in the context of our paper, might be true when τ is finite. Then we can expect to find a zero crossing in the close vicinity of the measurement time t (i.e., at a distance of the order of τ from t). However, when 0 < α < 1 the situation is more subtle. In this case the points on time axis will be clustered, and while their number is infinite visualizing these dots with a simulation we will observe a fractal dust.
In this Appendix we investigate the statistics of duration of the excursion straddling time t, deriving some known results along the way. Our approach is based on renewal theory, using methods given by Godreche and Luck [47]. The previous mathematical approaches are based on direct analysis of Brownian and Bessel motion, while we use the trick introduced in the main text. As mentioned, we replace the continuous regularized Bessel process or Brownian path, with a non-continuous one, with jumps of size after zero crossing, and at the end is taken to zero.
While α = (1 + D)/2D for the problem in the main manuscript, hence α ≥ 1/2, here we will assume that 0 < α < 1. The case α = 1/2 is the Brownian case, since a particle starting on will return for the first time to the origin according to the well-known law g(τ ) ∝ τ −3/2 . In our physical problem, the trick means that g * is dependent, Eq. (59). As in the main text, we denote {t 1 , t 2 , · · · , t k , · · ·} as times of the renewal events (i.e., in our problem zero crossing events). This is a finite set of times if is finite. The waiting times {τ 1 , τ 2 , · · ·} are independent identically distributed random variables, due to the Markovian property of the underlying paths, and t 1 = τ 1 , t k = k i=1 τ i . We call t the observation time, and the number of renewals in (0, t) is denoted by n, so by definition t n < t < t n+1 (this is obviously true as long as is finite). The time B = t − t n is called the backward recurrence time [47], the time F = t n+1 − t is the forward recurrence time, and ∆ = t n+1 − t n is the time interval straddling t. In our problem, the backward recurrence time is the duration of the meandering process which starts at t n and ends at t. For Brownian motion and Bessel processes, where t n and t n+1 are the zero hitting times straddling t, the statistics of F, B and ∆ are non-trivial.
We now obtain the PDF of ∆, denoted d t (∆), using the methods in [47]. From the constraint t n < t < t n+1 we have where δ[· · ·] is the Dirac δ function, and I(t n < t < t n+1 ) = 1 if the condition in the parenthesis is true, otherwise it is zero. We use the double Laplace transform Since the waiting times {τ i } are mutually independent, identically distributed random variables we find using t n = n i=1 τ i e −stn =ĝ n (s), whereĝ(s) = ∞ 0 e −sτ g(τ )dτ is the Laplace transform of g(τ ). A similar result holds for the other expressions in Eq. (184). It is then easy to find and summing the geometric series we get When g(τ ) = exp(−τ ) we find in the limit of large t d t (∆) → ∆ exp(−∆) so the minimum of this PDF is at ∆ = 0 and ∆ = ∞ which is expected. For our more interesting case,ĝ(u) ∼ 1 − g * |Γ(−α)|u α for small u, hence we get in the small u and s limit, with an arbitrary ratio between them This is a satisfying result implying that the solution does not depend on g * and hence it does not depend on . Thus the artificial cutoff, , which was introduced merely as a mathematical tool, does not alter our final formulas. Let∆ = ∆/t and denote byd(∆) its PDF in the long time limit. With a useful inversion formula, given in the Appendix of [47], we invert Eq. (188) to the time domain and find d ∆ = sin πα π(∆) 1+α 1 which exhibits a discontinuity of its derivative at∆ = 1. Thus the PDF of 0 < ∆ < ∞ has a cusp at t, which allows for the identification of the measurement time t from a histogram of ∆. For Brownian motion, α = 1/2 Instead of fixing t, we may draw t e > 0 from an exponential PDF, R exp(−Rt e ), in such a way that the mean of t e , t e = 1/R, is very large, so that the number of renewals in the time interval (0, 1/R) is large. Similar to the previous case, we define ∆ e = t n+1 − t n , which is called the duration of the interval straddling an independent exponential time, and t n < t e < t n+1 . This case was treated rigorously by Bertoin et al. [31], for the Bessel process. Using the renewal theory approach, we now obtain their main result on the statistics of ∆ e , with a few hand waving arguments. From the definition of the Laplace transform,f (s) = ∞ 0 f (t e ) exp(−st e )dt e , we see that for a function f (t e ) depending on a random variable t e , the later being exponentially distributed with mean equal to 1/R, the averaged function is the Laplace transform of f (t e ) evaluated at s = 1/R followed by multiplication with R, namely f (t e ) = Rf (R). Inserting s = R in Eq. (187) followed by multiplication with R, we get the Laplace transform of the PDF of ∆ e exp(−u∆ e ) = [(u + R) α − u α ]/R α , where we used the small R and u limit and as usual exp(−u∆ e ) = ∞ 0 e −u∆e d e (∆ e )d∆ e where d e (∆ e ) is the PDF of ∆ e . This is the known result for the Bessel process when scaled properly, i.e., R = 1 in [31] who also give the inverse Laplace transform of Eq. (191), thus providing an explicit formula for the PDF d e (∆ e ).
To see the connection between renewal theory, and statistics of the duration of Bessel excursion straddling an exponential time, notice that in [31] a Bessel process in dimension 0 < d < 2 is considered with the relation d = 2(1 − α). In our work we consider motion in a logarithmic potential in one dimension, which is easily mapped onto a Bessel process in dimension d = 1 − 1/D. Thus, the exponent α in Ref. [31] is the same as ours since as we have shown α = (1 + D)/(2D), Eq. (30). Notice that for optical lattices D > 0, hence this system is a physical example for a regularized Bessel process in dimension −∞ < d < 1. Finally, note that the original Bessel process considers the distance |r| from the origin of a Brownian motion in d dimension, this being non-negative |r| ≥ 0. Hence the Bessel process does not exhibit zero crossings, so the points on the time axis are zero hitting points, not zero crossing points. This is a minor technical issue, due the symmetry of the binding effective potential discussed in the manuscript (i.e., negative and positive excursions in the logarithmic potential are statistically identical).
Similarly, we can obtain the limiting distributions of B and F using renewal theory [47]. The PDF ofB = B/t (here t is fixed) is Inserting α = 1/2 we get the result forB, for zero crossing of Brownian motion, obtained by Chung (see Eq. 2.22 there). One can also quantify other aspects of zero crossings. For example, it is well known that the averaged number of renewals follows n ∼ t α /g * |Γ(−α)| and hence using Eq. (59) n ∝ t α / which diverges in the continuum limit → 0, as it should. The distribution of n/ n is the well known Mittag-Leffler distribution, with index α and unit mean.
To conclude we see that the backward and forward recurrence time scale linearly with t, and they exhibit non-trivial behavior, which can be obtained either from analysis of Brownian or Bessel processes, or using renewal theory with the trick. The latter is a very simple approach, which requires some basic results in renewal theory. Importantly, given the tools of renewal theory, the exponent α of the first passage time PDF g(τ ) ∼ τ −1−α which is investigated in Appendix A, determines uniquely the statistics of F, B and ∆, in the long measurement time limit and when α < 1. The fact that the backward recurrence time is long in the sense that it scales with the observation time t, explains why the mean square displacement x 2 , we have found in the main text, depends on the properties of the meander.