The velocity distribution of a classical gas of atoms in thermal equilibrium is the normal Maxwell distribution. It is well known that for sub-recoiled laser cooled atoms, Lévy statistics and deviations from usual ergodic behavior come into play. In a recent letter, we showed how tools from infinite ergodic theory describe the cool gas. Here, using the master equation, we derive the scaling function and the infinite invariant density of a stochastic model for the momentum of laser cooled atoms, recapitulating results obtained by Bertin and Bardou [Am. J. Phys. **76**, 630 (2008)] using life-time statistics. We focus on the case where the laser trapping is strong, namely, the rate of escape from the velocity trap is *R*(*v*) ∝ |*v*|^{α} for *v* → 0 and *α* > 1. We construct a machinery to investigate time averages of physical observables and their relation to ensemble averages. The time averages are given in terms of functionals of the individual stochastic paths, and here we use a generalization of Lévy walks to investigate the ergodic properties of the system. Exploring the energy of the system, we show that when *α* = 3, it exhibits a transition between phases where it is either an integrable or a non-integrable observable with respect to the infinite invariant measure. This transition corresponds to very different properties of the mean energy and to a discontinuous behavior of fluctuations. While the integrable phase is described by universal statistics and the Darling–Kac law, the more challenging case is the exploration of statistical properties of non-integrable observables. Since previous experimental work showed that both *α* = 2 and *α* = 4 are attainable, we believe that both phases could also be explored experimentally.

## I. INTRODUCTION

Laser cooled atoms and molecules are important for fundamental and practical applications.^{1–4} It is well known that Lévy statistics describes some of the unusual properties of cooling processes.^{5–9} For sub-recoil laser cooling, a special atomic trap in momentum space is engineered. The most efficient cooling is found when a mean trapping time, defined more precisely below, diverges.^{7} In this sense, the dynamics is time-scale-free. The fact that the characteristic time diverges implies that the processes involved are non-stationary. Furthermore, in the physics literature, they are sometimes called non-ergodic. As is well known, ergodicity is a fundamental aspect of statistical mechanics.

Ergodicity implies that time and ensemble averages of physical observables coincide. This is found when the measurement time is made long compared to the time-scale of the dynamics. However, in the context of sub-recoiled laser cooled atoms, this time diverges, and hence, no matter how long one measures, deviations from standard ergodic theory are prominent. Given that lasers replace heat baths in many cooling experiments, what are the ergodic properties of the laser cooled atoms? In other words, what replaces the usual ergodic statistical framework? While previous work thoroughly investigated the distribution of momentum,^{7} we highlight the role of the non-normalized quasi-steady state. Our goal is to show how tools of infinite ergodic theory describe the statistical properties of the ensemble and corresponding time averages of the laser cooled systems.

Infinite ergodic theory was investigated by mathematicians^{10–12} and, more recently, in physics.^{13–23} It has a deep relation with weak ergodicity breaking found in the context of glassy dynamics.^{24,25} The terminologies, which might seem at first conflicting, will be briefly discussed toward the end of this paper. Infinite ergodic theory deals with a peculiar non-normalized density, describing the long time limit; hence, it is sometimes called the infinite invariant density. Previous works in the field of sub-recoil laser cooling^{7,26} foresaw this quasi-steady state. Hence, we start with a recapitulation exposing the meaning of the infinite density using a master equation approach. We show how to use this tool to investigate the ensemble and time averages of physical observables and discuss the fluctuations. In particular, we investigate the energy of the system. Since the atoms are non-interacting, in a classical thermal setting, the energy of the atoms per degree of freedom is *k*_{B}*T*/2. In addition, this is obtained from the Maxwell velocity distribution, which is, of course, a perfectly normalized density. We will show, among other things, that the energy of sub-recoiled gas is obtained under certain conditions with a non-normalizable invariant density. A transition is exposed in the statistical properties of the energy when the fluorescence rate *R*(*v*) ∝ |*v*|^{α}, in the vicinity of zero velocity, is controlled. Since experimental work demonstrates the capability of a variation of *α* from *α* = 2 to *α* = 4 and at least, in principle, with pulse shaping and Raman cooling^{7,27,28} to other values of *α*, the rich phase diagram of ergodic properties we find seems to us within the reach of experimental investigation.

In this work, we are influenced by advances in the statistical theory of optical experiments, in particular single molecule tracking. In this field, the removal of the problem of ensemble averaging^{29} led to new insights into applications and limitations of standard ergodic theory, for example, in the context of diffusion of single molecules in the cell^{25,30} and the power law distributed sojourn times of dark and bright states of blinking quantum dots.^{31,32} Similarly, here, the time averages of a single trajectory of an atom/molecule in the process of laser cooling are studied theoretically, as fits this special issue. As mentioned, the omnipresent power law distributed sojourn times are also found in laser cooled gases, and they describe the times in the momentum trap in the vicinity of small velocities. These are responsible for the emergence of the infinite ergodicity framework.

A summary of our main results was published in a recent letter.^{33} Here, we focus on the derivations of key results, numerical demonstrations of convergence to the infinite density, its relation to the scaling solution, and later a detailed study of the integrable and non-integrable phases; the latter is more mathematically challenging, while the former is universal, as is described by the Darling–Kac theorem. We present the model and analyze the distribution of velocities with a master equation; this gives both a scaling solution and the infinite invariant density; see Sec. II. Section III is devoted to a discussion of both ensemble and time averages in generality. We then focus on the energy of the system, developing tools for analyzing the corresponding paths; see Sec. IV. In Sec. V, we explore the Darling–Kac phase, where the observable of choice is integrable with respect to the non-normalized state corresponding to *α* < 3. We then study the non-integrable phase 3 < *α* in Sec. VI. We end with open questions, perspective, and a summary.

## II. INFINITE DENSITY AND SCALING SOLUTION

Let *v* > 0 be the speed of the atom under the influence of sub-recoil laser cooling. The stochastic process for *v*(*t*), presented schematically in Fig. 1, is described by the following rules.^{7} At time *t* = 0, draw the speed *v*_{1} from the probability density function (PDF) *f*(*v*). Momentum is conserved until the atom experiences a jolt due to the interaction with the laser field. Hence, the speed of the atom will remain fixed for time $\tau \u03031$. The PDF of $\tau \u03031$ conditioned on *v*_{1} is exponential,

Here, *τ*(*v*) is the mean lifetime, which depends on the speed of the particle. After time $\tau \u03031$, we draw a new speed *v*_{2} from *f*(*v*). The process is, then, repeated. Namely, we now draw the second waiting time $\tau \u03032$ from the PDF in Eq. (1) with an updated lifetime *τ*(*v*_{2}). This process is, then, renewed. Given *f*(*v*) and *τ*(*v*), we are interested in the speed of the particle at time *t*. The event of the change of velocity is called below a collision or a jump. Here, we analyze the velocity distribution using a master equation approach. A different elegant approach to the problem was considered by Bertin and Bardou^{26} using the dynamics of the lifetimes.

*Master equation.* Let $\rho v,t$ be the PDF of *v* at time *t*. The repeated cooling process, as described above, leads to an evolution of $\rho v,t$, which is governed by a master equation

containing gain and loss terms, respectively. The independence of the parent distribution *f*(*v*′) from the previous velocity value *v* leads to a factorization of the transition rates from *v* to *v*′,

Using Eq. (3) and the normalization $\u222b0\u221ef(v\u2032)dv\u2032=1$ gives

which identifies *R*(*v*) as a jump rate, the rate of leaving the state with velocity *v*, which is the inverse of the lifetime, i.e., *R*(*v*) = 1/*τ*(*v*). With these assumptions, the master equation [Eq. (2)] simplifies to

An invariant density $\rho *v$, which zeros the time derivative on the left-hand side of the master equation, is obtained easily from Eq. (2), i.e.,

leading to $R(v\u2032)f(v)\rho *v\u2032=R(v)f(v\u2032)\rho *v$, which is fulfilled if, e.g., the *v*− dependence on both sides is identical, i.e., if $f(v)=R(v)\rho *vC$, or

and *C* is some constant. Equation (6) represents the detailed balance condition that holds for the model of laser cooled atoms.^{7} Below, we promote two scenarios for the invariant solution of Eq. (5). The first is well known and it is the case when $\rho *v$ is normalizable, and the second is the case when it is not. We will see below that even if the normalization integral diverges, which happens easily, then non-normalizable density $\rho *v$ of Eq. (7) is still meaningful and can be used for calculating certain time and ensemble averages.

The integral in Eq. (5) represents the process where atom’s state *v*′ is shifted due to the laser–atom interaction and the new velocity *v* is drawn from the PDF *f*(*v*). In Ref. 7, a uniform model for *f*(*v*) was investigated, and hence (unless stated otherwise), in our examples below, we choose *f*(*v*) = 1/*v*_{max} for 0 < *v* < *v*_{max}; otherwise, it is zero. In equilibrium, we have

and the normalization is $Z=\u222b0\u221e\tau (v)f(v)dv$. Here, the steady state exists in the usual sense and the normalized density *ρ*^{eq}(*v*) can be used to predict the ensemble and the corresponding time averages of the process. In particular, Birkhoff’s^{34} ergodic theory states that for an observable $O[v(t)]$, the time average is equal to the ensemble average,

where the ensemble average in equilibrium is

However, if

and *α* > 1, the above standard framework does not work since *Z* diverges. This is precisely the situation for laser cooled atoms where *α* = 2 or *α* = 4, 6 depending on the specific atom–light interaction process.^{7,26} Clearly, once *v* becomes small, then the lifetime is very long. This is the widely discussed mechanism of sub-recoil cooling; the atoms once slowed down will have a very long lifetime and, hence, remain in the cold state for a long time. The atoms, thus, pile up close to the zero velocity. This phase, i.e., *α* > 1, is the case where infinite ergodic theory plays a vital role, as we will show.

To solve the problem, we consider the long time limit, and then, for *v* ≠ 0, we have^{7,26,33}

This is a quasi-steady state in the sense that the numerator is proportional to the usual *ρ*^{eq}(*v*) [Eq. (8)]. It is valid for *v* beyond a small layer around *v* ≃ 0, which itself shrinks with time to zero; see below. It is easy to check that Eq. (11) is a valid solution by insertion into Eq. (5). Indeed, the time derivative on the left-hand side gives $(\xi \u0303\u22121)\rho (v,t)/t$, which is vanishing when *t* → *∞*, and the right-hand side gives zero just like for ordinary steady states. Here, the exponent $0<\xi \u0303<1$ is called the infinite density exponent, and it will be soon determined and similarly for the constant *b*.

Equation (11) is invalid for very small *v* because it diverges at *v* → 0 in such a way that it is non-integrable since *α* > 1. Following Refs. 7 and 26, we seek a scaling solution

which describes the inner region of slow atoms and, hence, the cooling effect. Here, $v\u221d1/t\gamma \u0303$, and hence, the velocity is small since *t* is large. Of course, the inner and outer solutions [Eqs. (11) and (12)] must match, and we will exploit this to find the dynamical exponents of this problem $\xi \u0303$ and $\gamma \u0303$. The latter is what we call the scaling exponent. Using Eq. (11), we have *ρ*(*v*, *t*) ∝ *v*^{−α} for *v* → 0, and this must match the inner solution; hence, we have

Note that *g*(*x*) is integrable, namely, it can be normalized when *α* > 1.

We insert Eq. (12) into Eq. (5) to find an equation for *g*(*x*). Some thought is required with respect to the upper limit of integration that stretches to infinite velocities on the right-hand side of Eq. (5). However, the scaling function does not describe large velocities; in fact, velocities of order *v* ∝ *v*_{max} are modeled by the quasi-steady state, Eq. (11), as mentioned. Let us, however, pretend that we do not know this and see how it comes out of the master equation. We recall that *t* is large and replace the upper limit of integration over the velocity by an upper cutoff *ht*^{β}. We call *β* the cutoff exponent. Using Eqs. (5) and (12), we then find

where *g*′(*x*) is the derivative of *g* with respect to *x*. Clearly, the natural scaled variable is $x=t\gamma \u0303v$. We now realize that only the small *v* behavior of *τ*(*v*) determines the properties of the scaling solution and, hence, *g*(*x*). This, as already pointed out, is because the particles pile up close to the zero velocity, and hence, only the small *v* limit matters. In contrast, note that we need the full structure of *τ*(*v*) to describe the outer solution [Eq. (11)]. Hence, now we replace *τ*(*v*) → *cv*^{−α} in Eq. (14), and after a change of variables, we find

Recall the large *x* behavior of *g*(*x*) [Eq. (13)], which gives *g*(*x*)*x*^{α} → const, and hence, we may find the long time limit of the integral in Eq. (15) and get

Here, we used $limt\u2192\u221ef(x/t\gamma \u0303)=f(0)$, which is a constant not equal zero; furthermore, *f*(0) is swallowed in $N$ with other constants and, hence, is not presented explicitly in Eq. (16). Equation (16) should be time independent; hence, we can find the exponents of the problem, $\gamma \u0303\alpha =1$ and *β* = 0. The latter is expected since the cutoff of velocity is of order *t*^{0}, i.e., *v*_{max}. The constant $N$ is eventually determined from normalization; see below.

We now use the scaling exponent

and from now on, $\gamma \u0303=\gamma $. This exponent is in the range 0 < *γ* < 1, and it describes the fat tail PDF of the waiting time within the momentum trap, $\varphi 1(\tau \u0303)\u221d\tau \u0303\u22121+\gamma $ [see the details below; this PDF is found via averaging the exponential PDF $q(\tau \u0303|v)$ from Eq. (1) over the random velocities]. Here, we see that the mean waiting time diverges, and hence, as mentioned in the Introduction, the process is scale-free. Generally, infinite ergodic theory is related to such processes, where the mean of the microscopic time-scale diverges.

The scaling function is, thus, determined by the following equation:

The substitution *y*(*x*) = *xg*(*x*) helps simplifying this equation even further, and eventually, the solution reads^{7}

With the help of the generalized Dawson integral $F(p,x)=exp(\u2212xp)\u222b0xexp(yp)dy$, which according to Ref. 35 can be expressed by the confluent hypergeometric function *M*(*a*, *b*, *z*), the scaling function can be written as $g(x)=N\gamma c\gamma xF(1\gamma ,xc\gamma )=N\gamma M(1,1+\gamma ,\u2212x1/\gamma c)$, which agrees with the result in Ref. 7. $N$ is determined from normalization $\u222b0\u221eg(x)dx=1$, and from formula 7.612(1) from Ref. 36, we find

As announced, the infinite density exponent $\xi \u0303$ can be found by matching the inner and outer solution. Specifically, the left part of the outer solution [Eq. (11)] is $\rho (v,t)\u221dv\u2212\alpha /t1\u2212\xi \u0303$ is matched with the right part of the inner solution. Equations (12) and (13) give *ρ*(*v*, *t*) ∝ *v*^{−α}/*t*^{1−γ}, resulting in

Given $N$ in Eq. (20), we can determine *b* in Eq. (11). Using integration by parts, Eq. (19)^{35} yields $g(x)\u223cNc/x1/\gamma $ for *x* → *∞*. Hence, from the scaling solution [Eq. (12)], for large *x* = *vt*^{γ}, we have

This is matched to the small *v* behavior of the infinite density solution [Eq. (11)], which using Eq. (10) gives

Thus, we find the constant *b*, given by $bf(0)=N$.

To summarize, we find that^{33}

The function $Iv(v)$ is non-normalizable since $Iv(v)\u221dv\u22121/\gamma $ for *v* → 0 and 1/*γ* = *α* > 1. This is why it is called an infinite density. The non-normalizability is hardly surprising since we take a perfectly normalized PDF *ρ*(*v*, *t*) and multiply it by *t*^{1−γ}; hence, the area under the product will diverge in the long time limit. Still, the function $Iv(v)$ can be used to calculate ensemble and time averages, as we will show below.

Figures 2 and 3 demonstrate the main results of this section using finite time simulations of the process. In Fig. 2, simulations are presented for the uniform velocity PDF *f*(*v*), and we use *v*_{max} = 1, while for the rate function, we use *c* = 1. Unless stated otherwise, these parameters will be used in Figs. 1–11. In Fig. 3, we show the case where the parent velocity distribution is exponential. While the scaling solution is not sensitive to the shape of *f*(*v*), the infinite density is, as demonstrated in Figs. 2 and 3.

The distribution for very small velocity is described in Eq. (12) with the scaling function *g*(*x*), Eq. (19). This scaling solution as a standalone is not a sufficient description of a cooled system for the following reason. From Eqs. (12) and (13), we have *ρ*(*v*, *t*) ∝ *v*^{−α} for large *v*, and for now, let us focus on the realistic choice *α* = 2. We find the awkward situation that the second moment of *v*, namely, the mean kinetic energy, would be infinite due to the fat tail of the solution, which is not what we expect from a cold system. Thus, the scaling solution *g*(*x*) describes the pile up of particles close to the zero velocity, but it also exhibits an unexpected heavy power law tail for large velocities. The infinite density describing the outer region cures this problem in the sense that it describes the large velocity cutoff; see Figs. 2 and 3. From here, we see that for calculating, for instance, the mean kinetic energy, we need the infinite density $Iv(v)$. The technical details of this calculation will be presented below. Furthermore, while the scaling function cures the non-normalizable trait of the infinite density, which stems from its small *v* behavior, we also find that the complement is true: the infinite density cures the unphysical divergence of the kinetic energy, found using the scaling solution, which is due to its unphysical large *v* properties. In short, we need both tools to describe the velocity distribution.

Finally, there is experimental evidence for the scaling solution. According to Eq. (12), the width of the velocity distribution is proportional to *t*^{−γ} and, hence, to *t*^{−1/α}.^{27,37} Experimentally, one may control the shape of the rate function *R*(*v*) = 1/*τ*(*v*) in the vicinity of zero, and engineered experiments provided the values *α* = 2 and *α* = 4. Note that *α* = 2 corresponds to an analytical expansion of *R*(*v*) ∝ *v*^{2} around its minimum at *v* = 0. Hence, the theory [Eq. (12)] predicts that the full width at half maximum (FWHM) of the velocity PDF decays like *t*^{−1/2} for *α* = 2 and as *t*^{−1/4} for *α* = 4. Both these behaviors were, indeed, observed experimentally up to times of order 20 milliseconds, which, given the lifetime of the atom, is considered very long in this field. This cooling experiment from 1995, using cesium, reports below 3 nK temperatures, and due to the non-Gaussian nature of the momentum packet, temperature is defined by the width of the momentum distribution and not by its variance. This reality is also consistent with quantum Monte Carlo simulations.^{7} In a later experiment, the momentum distribution of helium was recorded in full agreement with the statistical method.^{7,38} As far as we know, so far, the infinite density was not explored experimentally.

## III. ENSEMBLE AND TIME AVERAGES

We now consider the long time limit of averaged observables, with the latter denoted as $O[v(t)]$, so the observable is a function of the stochastic process *v*(*t*). In usual statistical physics, both the ensemble average and the time average of such observables are investigated, and we do the same here. Examples are the indicator function $O[v(t)]=I[va<v(t)<vb]$ and the kinetic energy of the particle $O[v(t)]=v2(t)=Ek(t)$ and *m*/2 = 1, with *m* being the mass of the atom. The indicator function equals unity if the condition *v*_{a} < *v*(*t*) < *v*_{b} is true; otherwise, it is zero, namely, it is an observable that switches at random times between 1 and 0; thus, this observable represents a dichotomous process; see Fig. 1. The kinetic energy *E*_{k}(*t*) is an observable that needs no introduction. The ensemble average is

and hence, using Eq. (24), we find

In this sense, the infinite invariant density $Iv(v)$ replaces the standard invariant density *ρ*^{eq}(*v*). The formula is valid, provided that the integral is finite, and such an observable $O(v)$ is called integrable with respect to the infinite density. Examples are the indicator function and the kinetic energy,

and

The former shows that the number of particles in the interval *v*_{a} < *v* < *v*_{b} is shrinking in time, provided that 0 < *v*_{a}, and the latter indicates that the energy of the system is decaying similarly, like ⟨*I*(*v*_{a} < *v*(*t*) < *v*_{b})⟩∝⟨*E*_{k}(*t*)⟩∝ 1/*t*^{1−γ}. Since $Iv(v)\u221dv\u2212\alpha $ for *v* → 0, we see that the second integral exists only when *α* < 3; hence, the kinetic energy is an integrable observable only if *α* < 3. If this condition does not hold, then the kinetic energy is called a non-integrable observable and, then, other rules apply; see below. Roughly speaking, integrable observables cure the non-normalizable trait of the infinite density found for *v* → 0, and importantly, these integrable observables are very basic to physical systems, for example, the energy.

We now consider the time average of an integrable observable

In standard ergodic theory, $limt\u2192\u221eO\u0304(t)=\u27e8O\u27e9eq$. In the current case, the dimensionless variable $\u03d2=\gamma O\u0304(t)/\u27e8O(t)\u27e9$ remains random in the long time limit, and as shown below, it satisfies the Darling–Kac theorem,^{10} provided that the observable is integrable. Here, the mean of ϒ in the long measurement time limit is unity, and importantly, in that limit, its PDF is time independent. First, let us consider the ensemble mean, namely, we consider an ensemble of paths and average over time and then over the ensemble

Here, we switched the order of the time and ensemble averages. Considering the long time limit and using Eq. (24),

where we used 0 < *γ* < 1. Hence, we conclude^{33} that

Thus, we established a relation between the time average and the ensemble average. The latter is obtained using the infinite density $Iv(v)$, and thus, this invariant density is not merely a tool for the calculation of the ensemble average, but rather it also gives information on the time average. More precisely, as discussed below, ϒ exhibits universal statistics independent of the observable, provided that it is an integrable observable. The 1/*γ* factor in Eq. (31) stems from the time averaging, and for example, if *α* = 1/*γ* = 2, we have $\u27e8O\u0304(t)\u27e9\u223c2\u27e8O(t)\u27e9$, or since both sides of this equation actually go to zero, a more refined statement is $limt\u2192\u221e\u27e8O\u0304(t)\u27e9/\u27e8O(t)\u27e9=2$.

## IV. KINETIC ENERGY

We consider the time average of the kinetic energy of the atom

in the limit of long times. While we focus on a particular observable, the theory presented below is actually rather general. As mentioned, depending on the value of *α*, the observable *v*^{2}(*t*′) can be either integrable with respect to the infinite density or not. Hence, here, we will develop a general theory, for the time averages, valid, in principle, whether the observable is integrable or not. To do so, we use tools from random walk theory,^{39,40} in particular certain types of Lévy walks.^{41–43}

We focus on the numerator of Eq. (32), which is a functional of the velocity path,

Since we have no potential energy, $S(t)$ is the action, which is increasing with time and $E\u0304k(t)=S(t)/t$. Since the speed of particle is a constant between collision events,

Here, *v*_{1} is the velocity drawn at the start of the process, *v*_{2} is the velocity drawn after the first collision, etc. The times $\tau \u0303i$ are the times between collision events. *t*_{B}(*t*) is the time elapsing between measurement time *t* and the last event in the sequence, also called the backward recurrence time.^{44,45} Here, we have the constraint

See schematics in Fig. 4. Finally, *N*(*t*) is the random number of collisions until time *t*.

Recall that at each collision event, we draw *v*_{i} from the parent PDF *f*(*v*) and then the waiting time $\tau \u0303i$ from Eq. (1) $q(\tau \u0303|vi)=R(vi)exp[\u2212\tau \u0303R(vi)]$, which is the conditional PDF of the waiting time, given a velocity *v*_{i}. Hence, the waiting times are not identically distributed unless the rate is a constant. Furthermore, in this section, we will assume that

and zero otherwise. Thus, we consider a uniform distribution of the speed, and furthermore, the rate is

which is the inverse of the mean decay time *τ*(*v*) in Eq. (10). Here, *c* is a constant with units of time multiplied by speed to the power of *α*.

We may write $si=(vi)2\tau \u0303i$ and, similarly, $sB(t)=(vN(t)+1)2tB(t)$ and, then,

At first, this appears to be a problem of the summation of *N* random variables, which is classical in many fields, in particular in the theory of random walks. However, here, *N*(*t*) is random and is determined by the sequence of waiting times, which, in turn, are correlated with the jump size, i.e., the increments of *s*_{i} of $S(t)$ [Eq. (38)] and also *τ*_{i}s in Eq. (35) are constrained by the measurement time *t*. In particular, when *α* > 1, the statistics is very different from normal.

In what sense are we dealing here with a generalized Lévy walk? The simplest form of the Lévy walk deals with the displacement of a particle −*∞* < *X*(*t*) < *∞*, which is given by $X(t)=\u2211iN(t)vi\tau i\u0303+vN+1tB(t)$. Here, the velocity *v*_{i} is, say, either +1 or −1 with probability half (unbiased random walk), while the travel time PDF is fat tailed. In our study, *S*(*t*) > 0 plays a role similar to *X*(*t*) of a biased Lévy walk. The main difference is the following. For standard Lévy walks, the distribution of the flight times $\tau \u0303$ is independent of the value of the velocity, i.e., the joint PDF of velocity and travel times is a product of the marginal PDFs. Here, we have a different situation as the velocity of the atom after the jump event is correlated with the random time of flight by a law determined with *R*(*v*) = 1/*τ*(*v*). However, as stated in the abstract, the action *S*(*t*) is performing a generalized type of Lévy walk, which is analyzed below.

### A. PDF of action increments and waiting times

The joint PDF of $s,\tau \u0303$, and *v* is

where *f*(*v*) and $q(\tau \u0303|v)$ are defined in Eqs. (1) and (36), respectively. Here, the delta function comes from the definition of one step, i.e., given a specific velocity and waiting time, the increment of the action is fixed. Subscript 3 indicates that we are dealing here with three random variables. We already mentioned that the action increments and the waiting times are correlated, and hence, to solve the problem, we need the joint PDF of s and $\tau \u0303$. This is obtained from integration of Eq. (39) over *v*, which gives

when $0\u2264s\u2264vmax2\tau \u0303$, and if this condition is not valid, the joint PDF is equal zero. In the context of random walk theory, such joint distributions are investigated in the context of *coupled* continuous time random walks^{41,43,46–49} since the increment s and waiting times $\tau \u0303$ are correlated.

Furthermore, integrating over *s*, we get the marginal PDF of the waiting times,

where *γ*(·, ·) is the lower incomplete gamma function. For large $\tau \u0303$, we get a power law decay

Hence, if *α* > 1, the mean waiting time diverges and const = *c*^{1/α}Γ(1 + 1/*α*)/(*αv*_{max}). As is well known,^{24,25} the divergence of the mean waiting time signals special ergodic properties since no matter how long we measure, we cannot exceed the mean time, and hence, ergodicity in its usual sense is broken. At short waiting times, we have $\varphi 1(\tau \u0303)\u223cvmax\alpha /[(1+\alpha )c]$.

The marginal PDF of the action increment is found integrating Eq. (39) over $\tau \u0303$ and *v*. We use the notation that the argument in the PDF defines the variable under study, and in this case, it is easy to show that

Then, for the example *α* = 2, we find *ϕ*_{1}(*s*) = exp(−*s*/*c*)/*c*, namely, an exponential decay. More generally, the mean of the action increments $\u27e8s\u27e9=\u222b0\u221es\varphi 1(s)ds$ is an important quantifier, and Eq. (43) gives

Clearly, the mean diverges when *α* > 3, which is critical for our discussion. Note the peculiarity of the case *α* = 2 as the mean action increment is independent of *v*_{max}.

We now examine the distribution of *s* in some detail. For *α* > 2, we change variables according to *z* = *v*^{α−2}*s*/*c* in Eq. (43) and, then, the lower limit of the integral, namely, *v* = 0 corresponds to *z* = 0. Note that when *α* < 2, the lower limit of integration over *z* transforms to infinity. We find

For large s, the lower incomplete gamma function is a constant equal to the corresponding Gamma function; hence, for large s, we find the power law tail *ϕ*_{1}(*s*) ∝ *s*^{−1−1/(α−2)}, and from here, we see again that when *α* > 3, the mean action increment diverges. For example, for the experimentally relevant case *α* = 4, we have

which gives $\varphi 1(s)\u223cc\pi /(4vmax)s\u22123/2$ when *s* → *∞*.

For *α* < 2, the integral [Eq. (43)] can be solved similarly and the solution can be expressed in terms of the upper incomplete gamma function,

Now, for large s, the distribution *ϕ*_{1}(*s*) has an exponential cutoff for all *α* with 1 ≤ *α* ≤ 2, as can be seen from the asymptotic expansion of the upper incomplete gamma function Γ(*a*, *x*) ∼ *x*^{a−1} exp(−*x*), resulting for 1 ≤ *α* < 2 in

and for *α* = 2, we already saw that $\varphi 1(s)=1cexp(\u2212sc)$. For small s, one gets a finite value $\varphi 1(s\u21920)=v0\alpha \u22122c(\alpha \u22121)$ for all *α* > 1. The case *α* = 1 is special in this respect because here, *ϕ*_{1}(*s*) is given by

which blows up at small s like *ϕ*_{1}(*s*) ∼ [−ln(*s*/(*v*_{max}*c*)) − *γ*_{1}]/(*v*_{max}*c*), where *γ*_{1} ≃ 0.577 is the Euler–Mascheroni constant. We see that *α* = 1, which marks the transition between phases with and without a finite mean waiting time, also exhibits a transition from a blow up at the origin of *ϕ*_{1}(*s*) to the case where the PDF of increments is a constant at short *s*.

To summarize, we see that the critical value *α* = 1 marks the transition from a finite mean waiting to an infinite mean, while *α* = 3 marks the transition from a finite mean action increment to a diverging one. The critical value *α* = 3 is specific to the observable under study, namely, the energy of the atom. Considering an observable like *v*(*t*)^{q} and *q* ≠ 2 would yield other values of the critical exponent. Still, energy is a basic concept for thermodynamics, so we focus on this particular observable. In addition, *α* = 2 marks a transition from a power law to an asymptotically exponentially decaying distribution of action increments. However, this does not play an essential role in our investigations of ergodic properties of the process. Importantly, the joint PDF of s and $\tau \u0303$ does not factorize, and hence, we need to treat the correlations between these variables, which is what is done in Subsection IV B.

### B. Governing Montroll–Weiss equations for the distribution of the total action

We will now investigate the formalism, giving the PDF of the action $S(t)$ at time *t*, which is denoted as $P(S,t)$. In Subsection IV C, we consider some of its long time properties. This density is, of course, normalized according to $\u222b0\u221eP(S,t)dS=1$. From the relation $E\u0304k(t)=S(t)/t$ and, then, from the distribution of $S(t)$, we can gain insights into the ergodic properties of $E\u0304k(t)$.

Let $QN(S,t)dt$ be the probability that the *N*th transition event takes place in the small interval (*t*, *t* + d*t*) and the value of the action $S(t)$ switched to $S$. It is given by the iteration rule

with $Q0(S,t)=\delta (S)\delta (t)$ being the initial condition. Equation (48) describes the basic property of the process: to arrive in $S$ at time *t* when the previous collision event took place at time $t\u2212\tau \u0303$, the previous value of the action was $S\u2212s$. Solving this equation is made possible with the convolution theorem of the Laplace transform. Let

be the double Laplace transform of $Qn(S,t)$, where $u\u2194S$ and *p* ↔ *t* are Laplace pairs. Then, from the convolution theorem, we have

where $\varphi \u03022(u,p)$ is the double Laplace transform of $\varphi 2(s,\tau \u0303)$. The PDF $P(S,t)$ is, in turn, given by

Here, we sum over the number of events *N* in the time interval (0, *t*) and take into consideration that the measurement time *t* is in an interval straddled by two collision events, and we, further, integrate over the backward recurrence time $tB\u2032$. Finally, the statistical weight function is

Here, the waiting time for the next jump is larger than $tB\u2032$ since by definition, the next transition takes place at times larger than *t*. In this equation, we average over the speed, i.e., integrate over *v*, draw the waiting time from an exponential PDF with a velocity dependent rate, and take into consideration the fact that the last increment of the action is given by $s=tB\u2032v2$, leading to the delta function.

Now, we consider the double Laplace transform of $P(S,t)$, which is denoted by $P\u0302(u,p)$. We again use the convolution theorem, and Eqs. (50) and (51) give, after summing a geometric series,^{33}

In the context of random walk theory, such a formula is called the Montroll–Weiss equation,^{50} though typically instead of a double Laplace transform one invokes a Laplace–Fourier transform.^{46,49} In general, to invert such an expression to the $S,t$ domain is hard; however, certain long time limits can be obtained. In particular, from the definition of the Laplace transform, the mean of the action in Laplace space is

Hence, to get $\u27e8S(t)\u27e9$, we need to deal with a single inverse Laplace transform from *p* back to *t*. Using Eq. (53), we find

Note that the double Laplace transform $\varphi \u03022(0,p)$ evaluated at *u* = 0 is merely the Laplace $\tau \u0303\u2192p$ transform of the waiting time PDF $\varphi 1(\tau \u0303)$, so $\varphi \u03022(0,p)=\varphi \u03021(p)$. This, of course, comes from the fact that by integration of the joint PDF $\varphi 2(s,\tau \u0303)$ over *s*, we obtain the marginal PDF $\varphi 1(\tau \u0303)$. The second expression on the right-hand side of Eq. (55) is a contribution to the total action from the backward recurrence time and, in the long time limit and for ordinary processes, with a finite mean waiting time, can be neglected. To investigate ergodicity, one needs to go beyond the mean and consider the full distribution of $S(t)$ or at least its variance.

### C. Energy is an integrable observable when 1 < *α* < 3

We now obtain $\u27e8E\u0304k(t)\u27e9$ using two approaches, the first is based on the Montroll–Weiss machinery Eq. (53) and the second exploits the infinite invariant density. When *α* < 1, we have a standard steady state since, then, the mean waiting time is finite. We now consider the case 1 < *α* < 3. According to Eq. (42), the mean waiting time diverges, while Eq. (44) yields a finite mean action increment ⟨*s*⟩. Furthermore, the observable *v*^{2}, representing the energy, is integrable with respect to the infinite density. This is because $Iv(v)\u221dv\u2212\alpha $ for *v* → 0, and hence, the integral $\u222b0\u221eIv(v)v2dv$ is finite when *α* < 3.

We are interested in the long time limit of $\u27e8S(t)\u27e9$, and in Appendix A, we derive a rather intuitive equation,

with ⟨*N*(*t*)⟩ being the mean number of collisions in the time interval (0, *t*). The mean ⟨*N*(*t*)⟩ ∝ *t*^{1/α} increases sub-linearly because the mean time interval between collision events diverges. This can be justified with a hand waving argument as follows. When the mean $\u27e8\tau \u0303\u27e9$ is finite, we expect from the law of large numbers that $\u27e8N(t)\u27e9\u223ct/\u27e8\tau \u0303\u27e9$. However, in our case, when *α* > 1, the denominator diverges and is replaced with the effective mean, namely, $\u27e8N(t)\u27e9\u221dt/\u222b0tt\u2032\varphi 1(t\u2032)dt\u2032\u221dt1/\alpha $. Using well-known formulas from renewal theory, briefly recapitulated in Appendix A, we find

This gives the mean of the time averaged kinetic energy of the particles (mass *m*/2 = 1) $\u27e8E\u0304k(t)\u27e9\u223c\u27e8S(t)\u27e9/t$.

#### 1. Infinite ergodic theory at work

Infinite ergodic theory makes the calculation of ⟨*E*_{k}(*t*)⟩ or $\u27e8E\u0304k(t)\u27e9$ easy in the sense that one needs only the knowledge of the invariant density $Iv(v)$. Equations (24), (36), and (37) give

which is clearly non-normalizable. Then, the ensemble average [Eq. (26)]

and Eq. (31) give $\u27e8E\u0304k(t)\u27e9=\u27e8Ek(t)\u27e9/\gamma $. On the other hand, $\u27e8E\u0304k(t)\u27e9=\u27e8S(t)\u27e9/t$, and as expected, Eq. (57) divided by *γ* gives the same results as in Eq. (59). The calculation using infinite ergodic theory is straightforward and is essentially similar to the averaging we perform with ordinary equilibrium calculations, and in that sense, it provides a tool more convenient compared to the Montroll–Weiss random walk approach. Of course, the latter yields insights since it connects the mean of the observable to the number of renewals ⟨*N*(*t*)⟩. Note that as *γ* → 1/3, Eq. (59) exhibits a blow up, and furthermore, lim_{γ→1/3}*v*^{3−1/γ} = 1, and in this sense, the energy is switching to a behavior that is *v*_{max} independent; more generally, it is independent of the parent velocity PDF *f*(*v*). This marks the transition to the phase where the energy is no longer an integrable observable, as discussed in Sec. VI.

## V. FLUCTUATIONS DESCRIBED BY THE DARLING–KAC THEOREM

We consider the fluctuations of the time average of the energy for 1 < *α* < 3. Since $Ek\u0304=S(t)/t$, we investigate the distribution of $S(t)$ for long times, namely, $P(S,t)$, using its double Laplace transform [Eq. (53)]. The details of the calculation are presented in Appendix B. It is the custom to consider normalized random variables with unit mean, namely, instead of $S(t)$, we investigate $\u03d2=S(t)/\u27e8S(t)\u27e9$. This, by definition, is also $\u03d2=E\u0304k(t)/\u27e8E\u0304k(t)\u27e9$, and using Eq. (31), $\u03d2=\gamma E\u0304k(t)/\u27e8Ek(t)\u27e9$, and importantly, the denominator can be obtained by a simple phase space average of *v*^{2}, employing the infinite density. Then, asymptotically, for large *t*, the PDF of ϒ is time invariant, and according to Appendix B, Eq. (B9) is given by the Mittag–Leffler law^{33}

Here, *l*_{γ,1}(*t*) is the one-sided Lévy stable distribution with index *γ* such that *l*_{γ,1}(*t*) ↔ exp(−*p*^{γ}) are Laplace pairs. For the energy observables, this formula is valid for 1/3 < *γ* < 1.

While not proven here, the Darling–Kac theorem states that this result is valid for any observable integrable with respect to the infinite measure, namely, the PDF of $\u03d2=O\u0304(t)/\u27e8O\u0304(t)\u27e9=\gamma O\u0304(t)/\u27e8O(t)\u27e9$ is asymptotically given by ML(ϒ). We demonstrated this universality with a few examples, choosing as observable both the kinetic energy and the indicator function; see Figs. 5–8. In the limit *γ* → 1, the PDF ML(ϒ) approaches a delta function, corresponding to Birkhoff’s ergodic theorem, while in the opposite limit *γ* → 0, the PDF ML(ϒ) is exponential with mean unity (see Figs. 5–8), which explores this trend. It should be recalled that the limit *γ* → 0 is valid here only for the indicator function and not for the energy since the former is integrable in the whole domain 0 < *γ* < 1 (if *v*_{a} > 0), while the latter is integrable only in the interval 1/3 < *γ* < 1. Mathematically, the Mittag–Leffler distribution holds in the long time limit for ϒ = *N*(*t*)/⟨*N*(*t*)⟩, where *N*(*t*) is the random number of collisions/renewals in the process until time *t*. Roughly speaking, when ⟨*s*⟩ is finite, we have $S(t)\u2243N(t)\u27e8s\u27e9$, and hence, the statistics of the time averaged kinetic energy and $S(t)$ are given by the same law as the statistics of *N*(*t*). What is remarkable is that a similar behavior holds for any integrable observable and that the mean of the observable is easy to compute with the infinite density.

It is easy to find the long time limit of the moments of the process using Eq. (B8) and, then, inverting to the time domain. For example, $\u27e8S\u0302(p)\u27e9\u223c\u27e8s\u27e9p\u2212\gamma \u22121/(T1)\gamma $ for small *p*, which gives Eq. (57). Similarly, the EB parameter characterizing the fluctuations of the time average^{52} in the long time limit is given by

and EB stands for ergodicity breaking. Thus, when *γ* = 1, the fluctuations vanish, since then, we enter the ergodic phase when the mean trapping time is finite. Recall that for the energy observable, this equation is valid for *γ* > 1/3 since we assumed that ⟨*s*⟩ is finite.

## VI. KINETIC ENERGY—THE NON-INTEGRABLE PHASE 3 < *α*

When considering observables such as the kinetic energy, we have three types of behaviors,

where we used Eqs. (12) and (26). This, as mentioned, corresponds to cases where the mean waiting time is finite (1 < *γ*); the mean time diverges, but the mean action increment ⟨*s*⟩ is finite (1/3 < *γ* < 1); and finally the case where both diverge (0 < *γ* < 1/3). In the first case, standard ergodic theory holds; in the second, the Darling–Kac theorem is valid for the energy observable; and finally, we have a phase where the kinetic energy is non-integrable with respect to the infinite density 0 < *γ* < 1/3, and this is the case that is now treated. Note that both *g*(*x*) and *ρ*^{eq}(*x*) are perfectly normalizable distributions, so the intermediate phase 1/3 < *γ* < 1 is in that sense unique.

We did not present the derivation of Eq. (62) for 0 < *γ* < 1/3 since the result is rather intuitive. It means that in this regime, the contribution to the kinetic energy comes from the slow atoms where the scaling solution is valid. Technically, we consider $\u222b0\u221ev2\rho (v,t)dv$ and divide the integration into two parts: in the first, the density *ρ*(*v*, *t*) in the small *v* inner region is the scaling solution, and in the second, the density is approximated by the outer solution, namely, the infinite density. Then, after integrating over the observable *v*^{2}, we can show that in the long time limit, the former is the leading term. For the case 1/3 < *γ* < 1, the observable is integrable with respect to the infinite measure and, then, the opposite situation is found. When *γ* = 1/3 or *γ* = 1, one finds logarithmic corrections not treated here.

Using Eq. (19), we find

Here, the kinetic energy is independent of the specific details of *R*(*v*) besides *c* and *γ* = 1/*α*, namely, the parameters controlling the behavior of the rate of escape from the trap at small *v*. This is very different for the case 1/3 < *γ*. Furthermore, according to Eq. (63), the velocity PDF *f*(*v*) is not at all influencing the long time dynamics of the mean kinetic energy, though we assume that this distribution has finite moments, and support on the zero velocity.

The integral *J*_{γ} in Eq. (63) can be evaluated analogously to the calculation of the normalization of *g*(*x*) in Eq. (20). The integral converges only for *γ* < 1/3, and one obtains

where we can explicitly see the divergence at *γ* = 1/3 and that lim_{γ→0}*J*_{γ} → 1/3. Inserting this into (63), for the ensemble averaged kinetic energy, we get

From here, we obtain the expectation of the time average immediately,

Equations (65) and (66) confirm that in this phase, in contrast to the phases with *γ* > 1/3, the average energy is independent of the details of the parent velocity distribution *f*(*v*); especially here, it is independent of *v*_{max}. In Appendix C, we derive Eq. (66) using action (54).

The asymptotic relations between the ensemble averaged kinetic energy $Ek(t)$ and its time average $Ek(t)\u0304$ in the various phases can, therefore, be summarized as^{33}

To obtain the fluctuations of the time averaged energy, namely, to calculate the EB parameter [Eq. (61)], further work is required. We need to evaluate the second moment of the action, which in Laplace space is $\u27e8S\u03022(p)\u27e9=\u22022P\u0302(u,p)\u2202u2|u=0$, and then, using Eq. (53),

Unfortunately, all four terms are contributing to the small *p* limit under investigation. This means that unlike the case 1/3 < *γ*, here, the effect of the last interval in the sequence, which we called the backward recurrence time, is dominating the statistics of the time average of the energy. The details of the calculations are presented in Appendix C. Asymptotically, for 0 < *γ* < 1/3, we find

This expression clearly differs from the EB parameter found in the Darling–Kac phase 1/3 < *γ*, Eq. (61). The latter is universal in the sense that it is valid for any integrable observable. In contrast, here, the fluctuations are specific to a non-integrable observable, namely, the energy, provided that 0 < *γ* < 1/3. Note that when *γ* → 1/3 from below and above, EB = 2Γ^{2}(4/3)/Γ(5/3) − 1 ≃ 0.766, namely, Eqs. (61) and (69) match, so the EB parameter is a continuous function of *γ*, while its derivative is not. In the range 0 < *γ* < 1/3, the EB parameter has a minimum, an effect that we cannot explain intuitively. In Fig. 9, we present numerical results for the EB parameter vs 0 < *γ* < 1, comparing it to the analytical theory. A clear transition is observed when the energy observable switches from an integral to a non-integrable observable, a transition found when *γ* = 1/3.

In the limit *γ* → 0, we find from Eq. (69) EB = 4/5. To understand this limit, we note that

where unlike the definition in Eq. (69) here we have the ensemble averaged energy not the time averages. In this limit, we have a stagnation effect in the sense that in the measurement time *t*, the system remains in a particular, though, random velocity state for practically the whole duration of the process, namely, $Ek\u0304=\u222b0tEk(t)dt/t=v2\u2061t/t=v2$. Note that we consider here the limit where *t* is made long, and only then, *γ* → 0. It is easy to find ⟨*E*_{k}⟩ = ⟨*v*^{2}⟩ and $\u27e8(Ek)2\u27e9=\u27e8v4\u27e9$ in the limit *γ* → 0 using Eqs. (19) and (20). In particular, exp(−*x*^{1/γ}/*c*) = 1 if 0 < *x* < 1; otherwise, it is zero since *γ* → 0. Hence, using Eq. (19), $\u27e8v2(t)\u27e9=\u222b01x\u222b0xdzdx(c/t)2\gamma =(1/3)(c/t)2\gamma $ and, similarly, ⟨*v*^{4}(*t*)⟩ = (1/5)(*c*/*t*)^{4γ}, and indeed, lim_{γ→0}EB = [(1/5) − (1/3)^{2}]/(1/3)^{2} = 4/5. To conclude, Eq. (69) makes perfect sense in the limit *γ* → 0, marking the stagnation of the dynamics, and *γ* → 1/3, marking the transition to the Darling–Kac phase; in between, the solution is not trivial.

## VII. DISTRIBUTION OF TIME AVERAGES IN THE NON-INTEGRABLE PHASE

We will now obtain the PDF of the random variable $\u03d2=S(t)/\u27e8S(t)\u27e9=E\u0304k(t)/\u27e8E\u0304k(t)\u27e9$, namely, the normalized time averaged kinetic energy, in the phase when this observable is non-integrable with respect to the infinite density. Recall that when the energy is integrable, we obtain the universal Mittag–Laffler law [Eq. (60)]. Unlike the latter case, the PDF of ϒ denoted as *P*_{α}(ϒ) will now depend on the microscopical details of the model, in particular the PDF of the speed after collisions *f*(*v*). Here, we find *P*_{α}(ϒ) for the model under study, namely, the case where *f*(*v*) is uniform. The analysis does not allow us to obtain a general solution for *P*_{α}(ϒ) for all 3 < *α*, and we mainly revert to approximations. This is unlike the variance, given in terms of the EB parameter, Eq. (69), which was calculated exactly. Since the calculations are lengthy, here, we provide the outline of the theory, focusing on three cases, *α* = 4, *α* = 6, and *α* → *∞*. Comparing the semi-analytical solution to simulations, we gain insight into a new type of transition, which shows up as a sudden blow up of *P*_{α}(ϒ) for ϒ → 0. The effect is certainly not found for the Mittag–Leffler distribution within the integrable phase *α* < 3.

The double Laplace transform of the action propagator *P*(*S*, *t*) is given by $P\u0302(u,p)$ and the Montroll–Weiss-type equation [Eq. (53)]. The technical problem is to invert this solution in the limit of long times corresponding to the Laplace variable *p* → 0 being small. The functions $P(S,t)$ and $P\u0302(u,p)$ attain scaling forms, denoted as $P(S,t)\u223ct1\u2212\alpha /2f\alpha (S/t1\u22122/\alpha )$ and *P*(*u*, *p*) ∼ (1/*p*)*g*_{α}(*u*/*p*^{1−2/α}). Here, the limit under study is *t* → *∞* and $S\u2192\u221e$, the ratio $S/t1\u22122/\alpha $ remaining fixed and similarly in Laplace space. Note that we showed in Eq. (C8) that $\u27e8S(t)\u27e9\u221dt1\u22122/\alpha $; hence, the scaling of $S$ with time we use here is consistent with that observation. The two scaling solutions are related by the laws of Laplace transform,

Substituting $x=St1\u22122/\alpha $ and setting *p* = 1 turn Eq. (71) into the integral equation

relating *g*_{α}(*y*) and *f*_{α}(*x*) via this convolution transform with the kernel

For the scaling function *g*_{α}(*y*), with *α* > 3, i.e., in the non-Darling–Kac phase, we obtain from the *p* → 0 limit of $P\u0302(u,p)$ [Eq. (53)] the following exact form:

The goal is to obtain from this by inversion of Eq. (72) the scaling function *f*_{α}(*x*), which is a rescaled version of the limiting probability density *P*_{α}(ϒ) of the normalized time average ϒ. Such an inversion can be achieved, in principle, by a Mellin transform of both sides of Eq. (72), resulting in (Ref. 64, pp. 997),

where

is the Mellin transform of *g*_{α}(*y*) and $K\alpha \u0303(s)$ and $f\u0303\alpha (s)$ is defined analogously. Solving Eq. (75) for $f\u0303\alpha (s)$ and applying the inverse Mellin transformation give

and finally, by rescaling, the desired limit density is obtained as

where

The rescaling implies that the mean is

as requested. While, in principle, along these steps the problem of calculating the limit distribution *P*_{α}(ϒ) is solved, a fully analytic solution is available only in two cases, *α* = 4 and *α* = *∞*. For *α* = 4, the integrals in Eq. (74) can be evaluated by the residue calculus, yielding after some calculations the simple result

with the Mellin transform $g\u03034(s)=\Gamma (s)\Gamma (1\u2212s)$. Since the Mellin transform $K\alpha \u0303(s)$ of the integral kernel is also known in full generality, $K\u0303\alpha (s)=\Gamma (s)\Gamma (1\u2212(1\u22122/\alpha )s)$, we get the quotient in Eq. (76), and in addition, we can invert from Mellin space to obtain $f4(x)=1\pi e\u2212x24$. The scaling factor $C4=xf4$ follows from the general relation between the *n*th derivative $g\alpha (n)(y=0)$ and the moments $xnf\alpha $ of *f*_{α}(*x*),

which follows directly from Eq. (72). For *α* = 4, we get $C4=xf4=2\pi $ so that for *P*_{4}(ϒ), according to Eq. (77), we obtain a half Gaussian distribution

as an exact result. It is merely a coincidence that this half Gaussian is also found in the Mittag–Leffler phase when *α* = 2; see Fig. 6. We can proceed similarly for *α* → *∞* because all integrals and Mellin transforms are also known exactly in this case, yielding

and eventually,

where we use the step function. This result diverges at the origin ϒ → 0, and we will soon show that this is valid for any *α* > 4. Note that PDF is cutoff sharply at ϒ = 3, an effect that is related to the underlying uniform distribution of *v*, Eq. (36). More precisely, we may explain this result by noting that the atom maintains a constant velocity for practically all the duration of the experiment since *α* → *∞* or *γ* = 1/*α* → 0. Then, ϒ = *v*^{2}/⟨*v*^{2}⟩, and using the uniform PDF of velocities Eq. (36), we get Eq. (82). In Fig. 10, we present simulation results for *α* = 50 and compare them to theory *α* → *∞*. These match well with the limiting PDF *P*_{∞}(ϒ), the exception is that the histogram is smeared and does not show the step-like structure of the limiting PDF, which is found at ϒ = 3.

For also the experimentally relevant case *α* = 6, we can still get from Eq. (74) by the residue calculus a fully analytic expression for *g*_{6}(*y*), but the result is very lengthy involving third roots, etc., and it does not simplify as in the case *α* = 4. Therefore, one cannot analytically calculate its Mellin transform. This led us to find an approximate scaling function $g6\u2248(y)$, which deviates only little from the true function *g*_{6}(*y*), but for which the Mellin transform is known analytically. The general idea is to match in $g6\u2248(y)$ the true small *y*-behavior and the large *y*-asymptotics of *g*_{6}(*y*), which we know analytically from an analysis of Eq. (74). The simple form

shares with the exact function *g*_{6}(*y*) the identical first and second derivative at *y* = 0 and the exponent of the asymptotic behavior for *y* → *∞*. The relative deviation of $g6\u2248(y)$ from *g*_{6}(*y*) is strictly less than 4.2%. This function $g6\u2248(y)$ can be Mellin transformed^{65} and leads via Eq. (76) and subsequent rescaling to an analytical expression for $P6\u2248(\u03d2)$, which can be expressed in terms of a Fox-H function^{66} as

### A. Accumulation effect for *α* > 4

We mentioned already that *P*_{α}(ϒ) diverges at the origin ϒ → 0 for any *α* > 4 and explained that this means that a large population of particles remain slow for the whole duration of the experiment. Now, we wish to characterize this effect more precisely.

It is easy to see that the exact asymptotic behavior of *g*_{α}(*y*) for *y* → *∞* is given by

with the decay exponent

and decay amplitude

From that, the small *x*-asymptotics of *f*_{α}(*x*) exactly follow as

with the exponent

and amplitude

The following table gives an overview for the resulting behavior of $f\alpha (x\u21920)\u223ct(\alpha )x\u2212s0(\alpha )$ and of the limit distribution *P*_{α}(ϒ → 0) = *C*_{α}*f*_{α}(*C*_{α}ϒ → 0). The results are obtained with the explicit form of *C*_{α} derivable from the exact form of $g\alpha \u2032(0)$,

This table shows that *α* = 4 takes a special role as it separates the diverging behavior of *f*_{α}(*x*) near *x* = 0 for *α* > 4 from the vanishing behavior for 3 < *α* < 4. Note also the diverging amplitude in *P*_{3+ɛ}(*x*) as *ɛ* → 0 in contrast to the vanishing amplitude of *f*_{3+ɛ}(*x*), which is a consequence of the diverging scaling factor *C*_{3+ɛ} ∼ *ɛ*^{−1}.

α
. | f_{α}(x → 0)
. | P_{α}(x → 0)
. |
---|---|---|

α → ∞ | $12x\u221212$ | $123x\u221212$ |

50 | $24\u2061sin(\pi 48)\pi \Gamma (2548)x\u22122348$ | $24\u2061sin(\pi 48)\pi \Gamma (2548)[(1+2\u2061cos(\pi 25))\Gamma (4925)]2548x\u22122348$ |

6 | $2\pi \Gamma (34)x\u221214$ | $1214\pi \Gamma (34)\Gamma (53)34x\u221214$ |

5 | $334\pi \Gamma (56)x\u221216$ | $33[12(3\u22125)]5124\pi \Gamma (56)\Gamma (85)56x\u221216$ |

4 | $1\pi x0$ | $2\pi x0$ |

$72$ | $338\pi \Gamma (76)x16$ | $338\pi \Gamma (76)cos(3\pi 14)sin(\pi 7)\Gamma (107)76x16$ |

$103$ | $23\pi \Gamma (54)x14$ | $(1+5)5212\pi \Gamma (54)\Gamma (75)54x14$ |

α = 3 + ɛ | $\epsilon x12\u2212\epsilon $ | $(23\pi \Gamma (75))32\epsilon \u221212x12\u2212\epsilon $ |

α
. | f_{α}(x → 0)
. | P_{α}(x → 0)
. |
---|---|---|

α → ∞ | $12x\u221212$ | $123x\u221212$ |

50 | $24\u2061sin(\pi 48)\pi \Gamma (2548)x\u22122348$ | $24\u2061sin(\pi 48)\pi \Gamma (2548)[(1+2\u2061cos(\pi 25))\Gamma (4925)]2548x\u22122348$ |

6 | $2\pi \Gamma (34)x\u221214$ | $1214\pi \Gamma (34)\Gamma (53)34x\u221214$ |

5 | $334\pi \Gamma (56)x\u221216$ | $33[12(3\u22125)]5124\pi \Gamma (56)\Gamma (85)56x\u221216$ |

4 | $1\pi x0$ | $2\pi x0$ |

$72$ | $338\pi \Gamma (76)x16$ | $338\pi \Gamma (76)cos(3\pi 14)sin(\pi 7)\Gamma (107)76x16$ |

$103$ | $23\pi \Gamma (54)x14$ | $(1+5)5212\pi \Gamma (54)\Gamma (75)54x14$ |

α = 3 + ɛ | $\epsilon x12\u2212\epsilon $ | $(23\pi \Gamma (75))32\epsilon \u221212x12\u2212\epsilon $ |

### B. Physical consequences

We discovered for *α* > 4 an accumulation effect, namely, the divergence of the PDF of the time averages, found at low energies, e.g., where *α* = 6 and ϒ → 0. This means that a significant population of atoms remains at small speeds for the whole duration of the experiment. In turn, this is useful when one wishes to reduce scattering or spatial spreading, namely, holding atoms close to the dark zero momentum state for long durations. Thus, while for the optimization of the commonly used relaxation time of the full width at half maximum (FWHM) of the velocity packet, which decays as *t*^{−1/α}, one should consider small values of *α* to obtain fast relaxation (say *α* = 2), to maintain some of the population with small kinetic energy for long durations, large values of *α* (say *α* = 6) are beneficial, as the trapping times become statistically longer. Surprisingly, *α* = 4 marks a quantitative transition of the low energy statistics, which we discovered from the analysis of the time averages.

## VIII. PERSPECTIVE

The rate of escape from the velocity trap, *R*(*v*) ∝ *v*^{α} for small *v*, implies that for laser cooled systems, the mean escape time is infinite when 1 < *α* or, equivalently, *γ* < 1. From the point of view of cooling, this is an advantage in the sense that typical speeds are low (nano-Kelvin regime). At the same time, this leads to the applicability of infinite ergodic theory, including the non-normalizable measure. The system shares some features that are similar to glassy dynamics, in particular the trap model.^{24} In that model, we have a density of states *ρ*(*E*) = exp[−*E*/*T*_{g}]/*T*_{g}, where *E* > 0 are the trap depths and *T*_{g} is a measure of disorder. The system is in contact with a heat bath at temperature *T*. We will not go into the details of the anomalous dynamics in this model; however, we also point out that here, we encounter a non-normalizable state. The partition function is $Z=\u222b0\u221e\rho (E)exp(E/T)dE$ (here, *k*_{B} = 1), and hence, it diverges when *T*_{g} > *T*. This low temperature glassy phase also corresponds to the case where the mean trapping time diverges and where one finds anomalous kinetics. Thus, for both the sub-recoiled system and the trap model, we find diverging mean trapping times and also the blow up of the normalization of the usual steady state. Bouchaud described such systems as exhibiting weak ergodicity breaking since one has exploration of phase space although time and ensemble averages differ.^{24,67} Mathematicians use the term infinite ergodic theory since they realize that the non-normalizable measure is the key ingredient of the theory. Furthermore, this non-normalizable state is related to the normalized distribution *ρ*(*v*, *t*); see Eq. (24), and it is approached from a broad class of initial conditions. Thus, some actually call the dynamics ergodic, i.e., the term infinite ergodic theory implies that the dynamics is ergodic, while in the physics literature, others describe it as non-ergodic. In short, one should distinguish between the operational definition of ergodicity, time and ensemble averages coincide, and the fact that in the long time limit, a unique density is approached, be it normalized or not. We conclude that weak ergodicity breaking and infinite ergodic theory are deeply related. The statistical theory applies not only to models in a non-equilibrium setting such as laser cooled atoms but also to systems with a canonical Boltzmann–Gibbs measure even if the latter is not normalized.^{19}

What are the consequences of laser cooling? Remarkably, using Eq. (62), we conclude that the most efficient cooling, in the sense of the fastest relaxation of the mean energy, is found for *γ* = 1/3. Thus, the transition in the ergodic properties of the system investigated here, which takes place when *α* = 3 or *γ* = 1/3, is physically connected to the optimal cooling of energy. This is not a coincidence, namely, at the transition point, the time dependence changes, although the fact that this point is optimal seems to us as merely good luck. In contrast, for the FWHM of the velocity packet,^{27} we do not have such an optimum. Instead, as mentioned, it decays like *t*^{−1/α}, favoring small values of *α* for faster relaxation.^{27} Thus, the classification of an observable as either integrable (energy, *γ* > 1/3) or non-integrable (energy *γ* < 1/3, FWHM) with respect to the infinite invariant measure is crucial, both mathematically and physically. We should note that the FWHM is not a dynamical observable in the sense that it cannot be obtained as a functional of a single particle path.

In the context of sub-recoil laser cooling, our work raised a few questions, and here, we point out possible extensions.

It is a challenge to see if quantum Monte Carlo simulations

^{7}can be used to numerically investigate the non-normalizable state and the time and ensemble averages.Experimentally finding the infinite density is another obvious challenge. Single atom experiments yield direct insights into the trajectories and the time averages.

^{68–71}In our examples, we used simple forms for

*f*(*v*) and*R*(*v*). It is important to realize that our main results are generally valid, such as the applicability of infinite ergodic theory, although it is clear that details do depend on the microscopical behavior of these functions. In this context, we have recently considered^{55}other models, including the case where the process is not renewed after each jolt. The main results are left unchanged.What happens to such a gas of atoms in the presence of some binding field, e.g., a harmonic trap? What is the pressure of the gas? What will happen when we add interactions? Will that drive the system to a true thermal state?

When we considered time averages, the measurement starts when the process is initiated (lower limit of the time integral is zero). Instead, one may prepare the system at time

*t*= 0, then wait until a time*t*_{a}, and only then perform a measurement, i.e., time average in the interval (*t*_{a},*t*_{a}+*t*). In this case, we expect that statistical properties of the system will depend on the aging time*t*_{a}. One can, then, wonder whether the infinite measure will also play an important role under these conditions. In this regard, we may be optimistic, see the modification of the Darling–Kac theorem to the aging regime, in the context of deterministic dynamics.^{63}Sisyphus cooling is also described in terms of Lévy processes,

^{6,56–61}and infinite*covariant*densities were studied in this context.^{14,62}Hence, the statistical framework of non-normalized states is, indeed, widely applicable.^{8}However, for Sisyphus cooling, the physics is orthogonal to the current one. The Sisyphus friction force is vanishing for large*v*, and thus, the non-normalizable trait comes from the high speed particles.^{14}Here, we have the opposite situation that the rate is anomalous for small*v*. Technically, this is related to the fact that for optical lattices, infinite covariant densities are studied, while here, the focus is on infinite invariant densities.According to the model, the width of the velocity distribution shrinks with time. In addition, as mentioned in the text, this was, indeed, observed in experiments until times of order milliseconds. Setting aside experiments, in the limit of

*t*→*∞*, the model predicts a complete pile up at the zero velocity, which seems to be a far-fetched idea. We have, thus, considered an idealized situation; in fact, one could introduce some cutoffs to the process using*R*(*v*) ∝*v*^{α}+*ϵ*, where*ϵ*is very small. Such a cutoff could be important as it would mean that in the very long time limit, the system will eventually relax to a normalized state. Indeed, in the absence of cutoffs, the time renewal process is scale-free, and hence, it is a random fractal. Like any fractal in nature, cutoffs could be important, as is well known.As shown here clearly and as well known, more generally, the theory of infinite ergodic theory is a theory of observables. For example, in our study, the indicator function

*I*(*v*_{a}<*v*(*t*) <*v*_{b}) is integrable with respect to the infinite density when 0 <*v*_{a}and non-integrable if*v*_{a}= 0. This is because of the non-integrable nature of the infinite density at small*v*, which, as stressed, is related to the fact that*R*(*v*) ∼*v*^{α}. The kinetic energy is integrable when 1/3 <*γ*= 1/*α*< 1, and this has consequences for the ergodic properties of the process. One could consider other observables such as |*v*|, and the main conclusions of this paper would be left unchanged. It should also be noted that when the invariant measure is finite and a usual steady state exists, an observable can be non-integrable, e.g., in a thermal setting, a particle in a harmonic trap has a Boltzmann density proportional to exp(−*kx*^{2}/2*k*_{B}*T*); hence, an observable that might seem a bit weird to some, such as $O(x)=exp(x4)$, is non-integrable. The case of non-integrable observables with respect to the steady state, i.e.,*α*< 1, is of some theoretical interest in the context of the stochastic model under study. At first, this might seem academic since so far, in the experiment,*α*> 1; however, this could be of interest in dimensions greater that unity; see below. An issue in our mind is whether a specific observable is*physically interesting or measurable*, and we worked under the assumption that energy is a physically worthy observable.We considered here the parent velocity distribution

*f*(*v*) being a constant at small*v*, e.g., a uniform velocity distribution. In Ref. 26, it was suggested that*f*(*v*) = const v^{d−1}for*v*<*v*_{max}and*d*is the dimension. Hence, as mentioned, the focus of this paper was on one-dimensional systems, mainly for the sake of simplicity.^{72}However, once again, the main conclusions of our paper are left unchanged, or more correctly, when minor adaptations are made, we may reach similar conclusions, for example, the infinite density in the general case, $Iv(v)\u221dv\u2212\alpha +d\u22121$, which is clearly non-normalizable when*α*−*d*+ 1 > 1. Hence, the cases*α*= 2 and*d*= 2 are special as it falls on the border between ergodicity in its usual sense and infinite ergodic theory. Such cases are left for future work.

## IX. SUMMARY

Our starting point was the master equation for the speed of the particle, which was previously studied with several methods.^{7,26} Here, we highlighted the infinite measure $Iv(v)$, which is a non-normalizable quasi-steady state of the system. To explore the ergodic properties of the system, we introduced a generalized Lévy walk approach. This tool gives a Montroll–Weiss-like formula, Eq. (53), which is a formal solution to the problem, but more importantly, it can be analyzed in the long time limit, giving statistical information on the distribution of the action $S(t)$ and, from it, the distribution of the time average of the energy $EK\u0304(t)=S(t)/t$. Here, we used the fact that between collision events, the momentum is conserved, so the speed is constant, changing abruptly at random times. With this method, we are able to obtain the properties of the time averages, which are functionals of the stochastic process. We focused on the kinetic energy of the particles; however, the approach presented here is more general. Technically, the increments of the walk are action increments s, and the joint PDF of the increments s and the waiting times are the basic ingredients of this coupled walk.

We find three phases in the ergodic properties of the process. The case *γ* > 1, corresponding to a finite mean time between collisions, was not considered here in detail since the standard ergodic theory holds as the invariant measure is normalizable. In the regime 1/3 < *γ* < 1, the Darling–Kac theorem holds for the observable of interest. This means that we may quantify the fluctuation of the time averages using the Mittag–Leffler law, Eq. (60), which is a universal type of statistical law in these types of problems. More precisely, this theorem is valid for observables that are integrable with respect to the infinite density. Finally, when 0 < *γ* < 1/3, the energy observable is non-integrable with respect to the infinite density. Here, the mean action increment ⟨*s*⟩ diverges together with the mean time between collisions. The consequence of these phases manifests themselves in several predictions. The decay with time of the ensemble energy, Eq. (62), goes through a qualitative change at the boundary between integrability and non-integrability *γ* = 1/3, similarly, for the relation between the mean of the time average and the ensemble average, Eq. (67). Finally, the *EB* parameter, Eqs. (61) and (69), characterizes the fluctuations of the time average, and it also exhibits a discontinuous behavior at *γ* = 1/3. Thus, we have exposed the rich consequences of the fact that an observable is tuned from being integrable to non-integrable. Interestingly, experiments use *α* = 4 (where energy is non-integrable) and *α* = 2 (where energy is integrable), so we believe that the classification we performed is of possible practical value. Finally, we discovered in Sec. VII another sort of transition. For *γ* < 1/4, the PDF *P*(ϒ) exhibits an accumulation effect, blowing up at ϒ → 0; see Fig. 11. This implies that some of the particles remain in the very cold phase, in the sense of very small velocities, for very long periods.

## ACKNOWLEDGMENTS

The support from the Israel Science Foundation under Grant No. 1614/21 is acknowledged (E.B.). This work was supported by the JSPS KAKENHI under Grant No. 240 18K03468 (T.A.). We thank Tony Albers, Nir Davidson, and Lev Khaykovich for helpful suggestions.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: RENEWAL THEORY

#### 1. Basics of renewal theory

We start with a brief recapitulation of renewal theory.^{44,45,73–75} Let $\varphi 1(\tau \u0303)$ be the PDF of time intervals between renewal events (in our case, collisions that modify the velocity). The process starts at time *t* = 0, and we draw a waiting time from the mentioned PDF, and this defines the point on the time axis for the first renewal event. We continue this way for the second event, etc. In the time interval (0, *t*), we have *N* events, and the latter is, of course, a random variable. Let *Q*_{N}(*t*)d*t* be the probability that the *N*th renewal event is taking place in the interval (*t*, *t* + *dt*). Then, from the renewal property of the process, we have

with the initial condition *Q*_{0}(*t*) = *δ*(*t*). The probability of finding *N* renewals in (0, *t*) is

Here, $W(t)=1\u2212\u222b0t\varphi 1(\tau \u0303)$ is the probability of not making a transition up to time *t*. Equation (A2), thus, describes a situation where the *N*th renewal takes place at time $t\u2212\tau \u0303$, and in the remaining time $\tau \u0303$, no jump was made. We now consider the Laplace transform of *P*_{N}(*t*) denoted as $P\u0302N(p)$, and using Eqs. (A1) and (A2) and the convolution theorem, we find

It follows that the mean of *N* is

To analyze the long time limit of ⟨*N*(*t*)⟩, we investigate $\u27e8N\u0302(p)\u27e9$ for small *p*. We are interested in the cases where the mean waiting time diverges, namely, following Eq. (41),

where in our case, $\tau 0=c[\Gamma (1+1/\alpha )/vmax]\alpha $. Then, in the limit *p* → 0, one can show that for *α* > 1,

where the leading term is the normalization. Inserting Eq. (A6) into Eq. (A4) and then performing a straightforward inverse Laplace transform, one finds $\u27e8N(t)\u27e9\u223c[\alpha \u2061sin(\pi /\alpha )/\pi ](t/\tau 0)1/\alpha $.^{44} This, in turn, gives the expression for ⟨*N*(*t*)⟩ in Eqs. (56) and (57) of the main text. The distribution of *N* in the long time limit is obtained using the small *p* behavior of Eq. (A3), and let us denote $\varphi \u03021(p)\u223c1\u2212b\gamma p\gamma $ and, then,

Inversion is made possible with the same tricks used to derive Eq. (B9),

where *l*_{γ,1}(·) is the one-sided Lévy PDF. Thus, the PDF of the action $S(t)$ discussed in the main text is the same as the PDF of the number of renewals *N* besides a scale, provided that 1/3 < *γ* < 1.

#### 2. Derivation of Eq. (56)

We now explain how to obtain Eq. (56) using Eq. (55). We investigate the latter in the limit *p* → 0 corresponding to long times. First, one can show that the second term on the right-hand side of Eq. (55) is negligible, provided that *α* < 3. Second, from the definition of the Laplace transform, we have

and using convolution, $\Phi \u03022(0,p)=[1\u2212\varphi \u03022(0,p)]/p$. Hence, from Eq. (55), we have

When *u* = 0, the Laplace transform of the joint PDF $\varphi \u03022(u,p)$ reduced to the Laplace transform of the marginal PDF of waiting times, namely, $\varphi \u03022(0,p)=\varphi \u03021(p)$. We now use Eq. (A4) and $\varphi \u03021(p)\u223c1$ to leading order in *p*, finding $\u27e8S\u0302(p)\u27e9\u223c\u27e8s\u27e9\u27e8N\u0302(p)\u27e9$, which gives Eq. (56).

### APPENDIX B: THE MITTAG-LEFFLER LAW

We consider the fluctuations of the time average of the energy for 1 < *α* < 3. Since $Ek\u0304=S(t)/t$, we investigate the distribution of $S(t)$ for long times, namely, $P(S,t)$, using its double Laplace transform [Eq. (53)]. From Eqs. (36), (37), and (39) and the definition of Laplace transform,

Integrating over s and $\tau \u0303$ gives

Similarly, we find

For *u* = 0, we have $\Phi \u03022(0,p)=[1\u2212\varphi \u03022(0,p)]/p$, which is the Laplace transform of the probability of not making a jump, also called the survival probability.

Since *t* is large and so is *S*, we consider the limits *u* → 0 and *p* → 0. This limit is considered under the condition that the ratio *p*^{γ}/*u* → const, or $S/t\gamma $, remains finite. This scaling is anticipated from the behavior of the moments of $S$, for example, Eq. (57).

We first consider the numerator of Eq. (53). Changing variables in Eq. (B3) according to *z*^{α} = *R*(*v*)/*p* = *v*^{α}/(*pc*), we find

As communicated, we consider the case 1/3 < *γ* < 1, and in this regime, *p*^{2γ−1}*u* = *p*^{3γ−1}(*u*/*p*^{γ}) → 0 since the ratio *u*/*p*^{γ} is fixed. We, then, take the upper limit of integration to infinity and, hence, to leading order

and $\u222b0\u221e(1+z\gamma )\u22121dz=\gamma \pi \u2061csc(\gamma \pi )$ with *γ* < 1.

We now consider the denominator of Eq. (53) and note that $limu,p\u21920\varphi \u03022(u,p)=1$ from normalization. Rewriting, we find

In the limit $G1\u223c(c\gamma p\gamma /vmax)\u222b0\u221edz(1+z\gamma )\u22121$ and $G2\u223cc(vmax)2\u2212\alpha u/(3\u2212\alpha )$ or using Eq. (44), *G*_{2} ∼ *u*⟨*s*⟩. We see that *G*_{1} ∝ *p*^{γ} and *G*_{2} ∝ *u*, so these two terms are of the same order.

with $(T1)\gamma =\gamma \pi c\gamma /(vmax\u2061sin\u2061\gamma \pi )$. Now, the inverse Laplace transform $u\u2192S$ gives

Let *l*_{γ,1}(*t*) be the one-sided Lévy stable distribution with index *γ* such that *l*_{γ,1}(*t*) ↔ exp(−*p*^{γ}) are Laplace pairs. Using −*d*/*dp* ↔ *t*, we obtain the PDF of the action for 1/3 < *γ* < 1, sometimes called the Mittag–Leffler law,

with $t\u0303=t/T1$ and $S\u0303=S/\u27e8s\u27e9$. The one-sided Lévy stable distribution is well documented,^{51} for example, in Mathematica, and hence, it is easy to plot this solution. More importantly, this distribution shows exactly how the fluctuations of the time averaged energy behave.

### APPENDIX C: MEAN ENERGY AND EB PARAMETER

#### 1. Formula for $\u27e8E\u0304k\u27e9$

We now derive Eq. (66). Specifically, we use as a model the uniform velocity PDF *f*(*v*) [Eq. (36)] and the rate function *R*(*v*) [Eq. (37)], although our final results are more general. The main tool is Eq. (55) for $\u27e8S\u0302(p)\u27e9$, which is evaluated in the limit of small *p* and then transformed to *t*, a standard Tauberian procedure valid in the long time limit. From here, as before, we get $\u27e8E\u0304k(t)\u27e9\u223c\u27e8S(t)\u27e9/t$.

Equation (55) is split into two terms,

with

When 1/3 < *γ*, the derivative $\u2202u\varphi \u03022(u,p)|u=0$ is the mean action increment ⟨*s*⟩ when *p* → 0. However, here, it diverges since we consider 0 < *γ* < 1/3. Recall that $\Phi \u03022(u,p)$ is associated with the probability of not jolting, namely, the term $H\u03022(p)$ is stemming from contribution to the total action from the backward recurrence time. In the limit, it was negligible for 1/3 < *γ*, but here, both terms are important.

To start, we rewrite Eq. (B2),

We set *u* to zero and change variables according to *v*^{α} = (*pc*)^{1/α}*z* and, then, get

The upper limit is then taken to be infinity, and we find in the small *p* limit $\varphi \u03022(0,p)\u223c1\u2212(pc)\gamma \pi \gamma /vmax\u2061sin\u2061\pi \gamma $. The first term here is the normalization, and the second indicates that the mean trapping time is diverging; more precisely, this well-known Tauberian relation comes from the fat tail of the waiting time PDF [Eq. (42)]. In these calculations and those that now follow, we use three integrals,

and

all valid in the regime of interest.

Using the same tricks, we find

and

Inserting Eqs. (C6) and (C8) into Eq. (C2), we find the small *p* limit of $H\u03021(p)$ and $H\u03022(p)$. It is now easy to transform from *p* to time *t* and find

The first term vanishes in the limit of *γ* → 0, namely, in that limit, the time interval (0, *t*) is effectively collision-free in such a way that the dominating contribution is the backward recurrence time^{53,54} [the *H*_{2}(*t*) term]. This means, roughly speaking, that in this limit, the backward recurrence time is equal to the measurement time. In contrast, when *γ* → 1/3, marking the transition to the integrability of the energy, namely, the Darling–Kac phase, the first term that stems from many jumps is diverging. From Eq. (C8), we get Eq. (66).

#### 2. EB parameter

To obtain the *EB* parameter, we use Eq. (68) to find the variance of $S$ in the long time limit. We consider the four terms *I*_{j}(*p*) with *j* = 1, …, 4 defined in Eq. (68) in the limit *p* → 0. The calculations are somewhat similar to those in Appendix C 1 although now we need to consider second order derivatives with respect to *u*; for example, we find

and using Eq. (C4),

Note that *I*_{1}(*p*) is independent of *v*_{max} in this limit, and so are the remaining terms *I*_{2}(*p*), *I*_{3}(*p*) and *I*_{4}(*p*), which are given by

and similarly,

Inverting to the time domain, we see that ⟨*S*^{2}(*t*)⟩∝ *t*^{2−4γ}. The integrals in Eqs. (C10)–(C13) are tabulated in Mathematica, so summing all the four terms in Eq. (68) and using ⟨*S*(*t*)⟩, Eq. (C8), we obtain the variance ⟨*S*^{2}(*t*)⟩ − ⟨*S*(*t*)⟩^{2}, and this after normalization yields the max *EB* parameter [Eq. (69)].

## REFERENCES

*et al.*, “

^{1}observables in one-dimensional intermittent maps with infinite invariant measures