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.

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 mathematicians10–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 cooling7,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 kBT/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 cooling7,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 averaging29 led to new insights into applications and limitations of standard ergodic theory, for example, in the context of diffusion of single molecules in the cell25,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.

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 v1 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 τ̃1. The PDF of τ̃1 conditioned on v1 is exponential,

(1)

Here, τ(v) is the mean lifetime, which depends on the speed of the particle. After time τ̃1, we draw a new speed v2 from f(v). The process is, then, repeated. Namely, we now draw the second waiting time τ̃2 from the PDF in Eq. (1) with an updated lifetime τ(v2). 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 Bardou26 using the dynamics of the lifetimes.

FIG. 1.

A sample path v(t) exhibits long sticking times whenever the particle is injected to a small velocity (a). This is due to the trapping mechanism, specifically the vanishing of the collision rate as R(v) = vα for small velocities. Here, we use α = 2 and the parent distribution f(v) is uniform in the interval (0, 1). We also show the indicator function that attains the value unity whenever 1/4 < v(t) < 1; otherwise, it is zero (b). As explained in the text, this observable is always integrable with respect to the infinite density, a trait crucial for infinite ergodic theory.

FIG. 1.

A sample path v(t) exhibits long sticking times whenever the particle is injected to a small velocity (a). This is due to the trapping mechanism, specifically the vanishing of the collision rate as R(v) = vα for small velocities. Here, we use α = 2 and the parent distribution f(v) is uniform in the interval (0, 1). We also show the indicator function that attains the value unity whenever 1/4 < v(t) < 1; otherwise, it is zero (b). As explained in the text, this observable is always integrable with respect to the infinite density, a trait crucial for infinite ergodic theory.

Close modal

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

(2)

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′,

(3)

Using Eq. (3) and the normalization 0f(v)dv=1 gives

(4)

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

(5)

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

(6)

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

(7)

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 ρ*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 ρ*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/vmax for 0 < v < vmax; otherwise, it is zero. In equilibrium, we have

(8)

and the normalization is Z=0τ(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’s34 ergodic theory states that for an observable O[v(t)], the time average is equal to the ensemble average,

(9)

where the ensemble average in equilibrium is

However, if

(10)

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 have7,26,33

(11)

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 (ξ̃1)ρ(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<ξ̃<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

(12)

which describes the inner region of slow atoms and, hence, the cooling effect. Here, v1/tγ̃, 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 ξ̃ and γ̃. 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

(13)

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 vvmax 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

(14)

where g′(x) is the derivative of g with respect to x. Clearly, the natural scaled variable is x=tγ̃v. 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

(15)

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

(16)

Here, we used limtf(x/tγ̃)=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, γ̃α=1 and β = 0. The latter is expected since the cutoff of velocity is of order t0, i.e., vmax. The constant N is eventually determined from normalization; see below.

We now use the scaling exponent

(17)

and from now on, γ̃=γ. This exponent is in the range 0 < γ < 1, and it describes the fat tail PDF of the waiting time within the momentum trap, ϕ1(τ̃)τ̃1+γ [see the details below; this PDF is found via averaging the exponential PDF q(τ̃|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:

(18)

The substitution y(x) = xg(x) helps simplifying this equation even further, and eventually, the solution reads7 

(19)

With the help of the generalized Dawson integral F(p,x)=exp(xp)0xexp(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γcγxF(1γ,xcγ)=NγM(1,1+γ,x1/γc), which agrees with the result in Ref. 7. N is determined from normalization 0g(x)dx=1, and from formula 7.612(1) from Ref. 36, we find

(20)

As announced, the infinite density exponent ξ̃ can be found by matching the inner and outer solution. Specifically, the left part of the outer solution [Eq. (11)] is ρ(v,t)vα/t1ξ̃ is matched with the right part of the inner solution. Equations (12) and (13) give ρ(v, t) ∝ vα/t1−γ, resulting in

(21)

Given N in Eq. (20), we can determine b in Eq. (11). Using integration by parts, Eq. (19)35 yields g(x)Nc/x1/γ for x. Hence, from the scaling solution [Eq. (12)], for large x = vtγ, we have

(22)

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

(23)

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

To summarize, we find that33 

(24)

The function Iv(v) is non-normalizable since Iv(v)v1/γ 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 t1−γ; 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 vmax = 1, while for the rate function, we use c = 1. Unless stated otherwise, these parameters will be used in Figs. 111. 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.

FIG. 2.

The scaled PDF of velocity t1−γρ(v, t) vs v for increasing times. A cooling process with a uniform parent velocity PDF f(v) and the maximum allowed velocity being unity is considered. The collision rate is R(v) = v2, so α = 2 and, hence, γ = 1/2. The data in the long time limit converge to the infinite density [Eq. (24)] presented with a black dashed line.

FIG. 2.

The scaled PDF of velocity t1−γρ(v, t) vs v for increasing times. A cooling process with a uniform parent velocity PDF f(v) and the maximum allowed velocity being unity is considered. The collision rate is R(v) = v2, so α = 2 and, hence, γ = 1/2. The data in the long time limit converge to the infinite density [Eq. (24)] presented with a black dashed line.

Close modal
FIG. 3.

Similar to Fig. 2, we show the scaled PDF of velocity t1−γρ(v, t) vs v; however, here, we simulated a process with an exponential parent velocity PDF f(v) with mean unity; the collision rate is unchanged, and hence, α = 2. The data in the long time limit converge to the infinite density [Eq. (24)], which is shown as a black dashed line. For finite times, we see expected deviations between the infinite invariant density and the numerical solution, which are obvious and clearly visible for small velocities. For not too large velocities, i.e., velocities much smaller than the mean of f(v), which is 1/2 in this example, the scaling solution, Eq. (12) with Eq. (19) (red line), nicely matches the data. We see that both the scaling solution and the infinite density describe the distribution of velocities for t = 104. For intermediate velocities, the two solutions match. Comparing with Fig. 2, we see that the sharp cutoff of the infinite density (due to the uniform PDF of velocities) is replaced here with an exponential decay.

FIG. 3.

Similar to Fig. 2, we show the scaled PDF of velocity t1−γρ(v, t) vs v; however, here, we simulated a process with an exponential parent velocity PDF f(v) with mean unity; the collision rate is unchanged, and hence, α = 2. The data in the long time limit converge to the infinite density [Eq. (24)], which is shown as a black dashed line. For finite times, we see expected deviations between the infinite invariant density and the numerical solution, which are obvious and clearly visible for small velocities. For not too large velocities, i.e., velocities much smaller than the mean of f(v), which is 1/2 in this example, the scaling solution, Eq. (12) with Eq. (19) (red line), nicely matches the data. We see that both the scaling solution and the infinite density describe the distribution of velocities for t = 104. For intermediate velocities, the two solutions match. Comparing with Fig. 2, we see that the sharp cutoff of the infinite density (due to the uniform PDF of velocities) is replaced here with an exponential decay.

Close modal

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) ∝ v2 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.

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 va < v(t) < vb 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 Ek(t) is an observable that needs no introduction. The ensemble average is

(25)

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

(26)

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

(27)

The former shows that the number of particles in the interval va < v < vb is shrinking in time, provided that 0 < va, and the latter indicates that the energy of the system is decaying similarly, like ⟨I(va < v(t) < vb)⟩∝⟨Ek(t)⟩∝ 1/t1−γ. Since Iv(v)vα 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

(28)

In standard ergodic theory, limtŌ(t)=Oeq. In the current case, the dimensionless variable ϒ=γŌ(t)/O(t) 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

(29)

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

(30)

where we used 0 < γ < 1. Hence, we conclude33 that

(31)

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 Ō(t)2O(t), or since both sides of this equation actually go to zero, a more refined statement is limtŌ(t)/O(t)=2.

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

(32)

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 v2(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,

(33)

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

(34)

Here, v1 is the velocity drawn at the start of the process, v2 is the velocity drawn after the first collision, etc. The times τ̃i are the times between collision events. tB(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

(35)

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

FIG. 4.

Action-time diagram illustrating the generalized Lévy walk performed by the action S(t), Eq. (38). tN(t) and SN(t) are the partial sums appearing in Eqs. (35) and (38), where tB(t) = ttN(t) and sB(t) = S(t) − SN(t) are the backward recurrence time and the associated action increment, respectively. For this sample trajectory, the number N(t) of renewals until time t is N(t) = 4.

FIG. 4.

Action-time diagram illustrating the generalized Lévy walk performed by the action S(t), Eq. (38). tN(t) and SN(t) are the partial sums appearing in Eqs. (35) and (38), where tB(t) = ttN(t) and sB(t) = S(t) − SN(t) are the backward recurrence time and the associated action increment, respectively. For this sample trajectory, the number N(t) of renewals until time t is N(t) = 4.

Close modal

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

(36)

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

(37)

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τ̃i and, similarly, sB(t)=(vN(t)+1)2tB(t) and, then,

(38)

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 si of S(t) [Eq. (38)] and also τis in Eq. (35) are constrained by the measurement time t. In particular, when α > 1, the statistics is very different from normal.

Remark.

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)=iN(t)viτĩ+vN+1tB(t). Here, the velocity vi 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 τ̃ 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.

The joint PDF of s,τ̃, and v is

(39)

where f(v) and q(τ̃|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 τ̃. This is obtained from integration of Eq. (39) over v, which gives

(40)

when 0svmax2τ̃, 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 walks41,43,46–49 since the increment s and waiting times τ̃ are correlated.

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

(41)

where γ(·, ·) is the lower incomplete gamma function. For large τ̃, we get a power law decay

(42)

Hence, if α > 1, the mean waiting time diverges and const = c1/αΓ(1 + 1/α)/(αvmax). 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 ϕ1(τ̃)vmaxα/[(1+α)c].

The marginal PDF of the action increment is found integrating Eq. (39) over τ̃ 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

(43)

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 s=0sϕ1(s)ds is an important quantifier, and Eq. (43) gives

(44)

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 vmax.

We now examine the distribution of s in some detail. For α > 2, we change variables according to z = vα−2s/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

(45)

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

(46)

which gives ϕ1(s)cπ/(4vmax)s3/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,

(47)

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) ∼ xa−1 exp(−x), resulting for 1 ≤ α < 2 in

and for α = 2, we already saw that ϕ1(s)=1cexp(sc). For small s, one gets a finite value ϕ1(s0)=v0α2c(α1) 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/(vmaxc)) − γ1]/(vmaxc), 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 τ̃ does not factorize, and hence, we need to treat the correlations between these variables, which is what is done in Subsection IV B.

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 0P(S,t)dS=1. From the relation Ēk(t)=S(t)/t and, then, from the distribution of S(t), we can gain insights into the ergodic properties of Ēk(t).

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

(48)

with Q0(S,t)=δ(S)δ(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τ̃, the previous value of the action was Ss. Solving this equation is made possible with the convolution theorem of the Laplace transform. Let

(49)

be the double Laplace transform of Qn(S,t), where uS and pt are Laplace pairs. Then, from the convolution theorem, we have

(50)

where ϕ̂2(u,p) is the double Laplace transform of ϕ2(s,τ̃). The PDF P(S,t) is, in turn, given by

(51)

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. Finally, the statistical weight function is

(52)

Here, the waiting time for the next jump is larger than tB 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=tBv2, leading to the delta function.

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

(53)

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

(54)

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

(55)

Note that the double Laplace transform ϕ̂2(0,p) evaluated at u = 0 is merely the Laplace τ̃p transform of the waiting time PDF ϕ1(τ̃), so ϕ̂2(0,p)=ϕ̂1(p). This, of course, comes from the fact that by integration of the joint PDF ϕ2(s,τ̃) over s, we obtain the marginal PDF ϕ1(τ̃). 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.

We now obtain Ēk(t) 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 v2, representing the energy, is integrable with respect to the infinite density. This is because Iv(v)vα for v → 0, and hence, the integral 0Iv(v)v2dv is finite when α < 3.

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

(56)

with ⟨N(t)⟩ being the mean number of collisions in the time interval (0, t). The mean ⟨N(t)⟩ ∝ t1/α 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 τ̃ is finite, we expect from the law of large numbers that N(t)t/τ̃. However, in our case, when α > 1, the denominator diverges and is replaced with the effective mean, namely, N(t)t/0ttϕ1(t)dtt1/α. Using well-known formulas from renewal theory, briefly recapitulated in  Appendix A, we find

(57)

This gives the mean of the time averaged kinetic energy of the particles (mass m/2 = 1) Ēk(t)S(t)/t.

1. Infinite ergodic theory at work

Infinite ergodic theory makes the calculation of ⟨Ek(t)⟩ or Ēk(t) easy in the sense that one needs only the knowledge of the invariant density Iv(v). Equations (24), (36), and (37) give

(58)

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

(59)

and Eq. (31) give Ēk(t)=Ek(t)/γ. On the other hand, Ēk(t)=S(t)/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/3v3−1/γ = 1, and in this sense, the energy is switching to a behavior that is vmax 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.

We consider the fluctuations of the time average of the energy for 1 < α < 3. Since Ek̄=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 ϒ=S(t)/S(t). This, by definition, is also ϒ=Ēk(t)/Ēk(t), and using Eq. (31), ϒ=γĒk(t)/Ek(t), and importantly, the denominator can be obtained by a simple phase space average of v2, 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 law33 

(60)

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 ϒ=Ō(t)/Ō(t)=γŌ(t)/O(t) 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. 58. 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. 58), 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 va > 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)N(t)s, 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.

FIG. 5.

To investigate the ergodic properties of the system, we simulate the cooling process model with a uniform velocity distribution f(v) with maximum velocity unity vmax = 1 and the rate R = vα, so c = 1. We, then, follow velocity paths as a function of time, and from each path, we obtain the time the particle spent in an interval va < v < vb. Namely, our observable is the indicator function. For standard ergodic processes, e.g., for the velocity of a gas particle obeying Maxwell statistics, this time divided by a long measurement time is approaching the probability of being in the mentioned velocity domain. For the laser cooling model, this time average fluctuates and the Mittag–Leffler distribution describes the corresponding statistics. Specifically, we show the histogram of the normalized random variable ϒ=0tI(va<v(t)<vb)dt/0tI(va<v(t)<vb)dt, which perfectly matches the theory [Eq. (60)]. Here, α = 1.25, so γ = 0.8, namely, we are not too far from the ergodic phase (γ = 1); hence, we see a peak in the distribution of ϒ close to its mean, which is unity. For an ergodic process, we will find a delta peak centered on unity, and this can be found in our model by choosing α < 1. We used t = 109, va = 0.4, and vb = 1.

FIG. 5.

To investigate the ergodic properties of the system, we simulate the cooling process model with a uniform velocity distribution f(v) with maximum velocity unity vmax = 1 and the rate R = vα, so c = 1. We, then, follow velocity paths as a function of time, and from each path, we obtain the time the particle spent in an interval va < v < vb. Namely, our observable is the indicator function. For standard ergodic processes, e.g., for the velocity of a gas particle obeying Maxwell statistics, this time divided by a long measurement time is approaching the probability of being in the mentioned velocity domain. For the laser cooling model, this time average fluctuates and the Mittag–Leffler distribution describes the corresponding statistics. Specifically, we show the histogram of the normalized random variable ϒ=0tI(va<v(t)<vb)dt/0tI(va<v(t)<vb)dt, which perfectly matches the theory [Eq. (60)]. Here, α = 1.25, so γ = 0.8, namely, we are not too far from the ergodic phase (γ = 1); hence, we see a peak in the distribution of ϒ close to its mean, which is unity. For an ergodic process, we will find a delta peak centered on unity, and this can be found in our model by choosing α < 1. We used t = 109, va = 0.4, and vb = 1.

Close modal
FIG. 6.

Same as Fig. 5; however, now, the observable is the kinetic energy v2(t) and α = 2. For this choice of α, the kinetic energy is integrable with respect to the infinite density, and hence, fluctuations of its time average obey the Darling–Kac law. We present the ML(ϒ) PDF of the normalized time average ϒ=0tv2(t)dt/0tv2(t)dt, Eq. (60). Since γ = 1/2, the Mittag–Leffler PDF is half a Gaussian presented as the theory here. We used 106 particles, and the measurement time was t = 109.

FIG. 6.

Same as Fig. 5; however, now, the observable is the kinetic energy v2(t) and α = 2. For this choice of α, the kinetic energy is integrable with respect to the infinite density, and hence, fluctuations of its time average obey the Darling–Kac law. We present the ML(ϒ) PDF of the normalized time average ϒ=0tv2(t)dt/0tv2(t)dt, Eq. (60). Since γ = 1/2, the Mittag–Leffler PDF is half a Gaussian presented as the theory here. We used 106 particles, and the measurement time was t = 109.

Close modal
FIG. 7.

Distribution of ϒ where the latter is the time spent by a path in a domain (va, vb) over its mean (the same as in Fig. 5); however, now, α = 3. Like any other observable integrable with respect to the infinite density, this normalized functional of the path obeys Mittag–Leffler statistics [Eq. (60)]. Since γ = 1/3, the peak of the histogram of ϒ is at the origin unlike the case presented in Fig. 5. In that sense, the fluctuations of time averages are larger when γ is small compared to γ close to unity as expected. When γ → 0, the distribution of ϒ for this integrable observable is exponential and, then, the corresponding EB parameter attains its largest value EB = 1.

FIG. 7.

Distribution of ϒ where the latter is the time spent by a path in a domain (va, vb) over its mean (the same as in Fig. 5); however, now, α = 3. Like any other observable integrable with respect to the infinite density, this normalized functional of the path obeys Mittag–Leffler statistics [Eq. (60)]. Since γ = 1/3, the peak of the histogram of ϒ is at the origin unlike the case presented in Fig. 5. In that sense, the fluctuations of time averages are larger when γ is small compared to γ close to unity as expected. When γ → 0, the distribution of ϒ for this integrable observable is exponential and, then, the corresponding EB parameter attains its largest value EB = 1.

Close modal
FIG. 8.

We show the dimensionless ratio Ō/O obtained from numerical simulations vs time. In the long time limit and if the observable is integrable with respect to the infinite density, we have, according to Eq. (31), Ō/O=α; a behavior we illustrate here for the indicator function, which as mentioned in the text, is an integrable observable for any choice of α. The energy observable is non-integrable for α = 4. Then, theory [Eq. (67)] predicts Ō/O=1/(12/α)=2, nicely matching the numerics.

FIG. 8.

We show the dimensionless ratio Ō/O obtained from numerical simulations vs time. In the long time limit and if the observable is integrable with respect to the infinite density, we have, according to Eq. (31), Ō/O=α; a behavior we illustrate here for the indicator function, which as mentioned in the text, is an integrable observable for any choice of α. The energy observable is non-integrable for α = 4. Then, theory [Eq. (67)] predicts Ō/O=1/(12/α)=2, nicely matching the numerics.

Close modal

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, Ŝ(p)spγ1/(T1)γ for small p, which gives Eq. (57). Similarly, the EB parameter characterizing the fluctuations of the time average52 in the long time limit is given by

(61)

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.

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

(62)

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 0v2ρ(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 v2, 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

(63)

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

(64)

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

(65)

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

(66)

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 vmax. 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)̄ in the various phases can, therefore, be summarized as33 

(67)

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 Ŝ2(p)=2P̂(u,p)u2|u=0, and then, using Eq. (53),

(68)

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

(69)

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.

FIG. 9.

The EB parameter vs γ for the kinetic energy observable v2. When EB = 0, we have an ergodic phase found when γ > 1. This ergodicity breaking parameter is a measure for the fluctuations of the time averages, and for 1/3 < γ < 1, it is given by Eq. (61). That formula is valid whenever the observable is integrable with respect to the infinite measure. For 0 < γ < 1/3, the kinetic energy observable is non-integrable, and the system enters a different phase, described in Eq. (69). The dots are obtained from finite time simulations of a model with a parent velocity distribution f(v), which is uniform, and the number of particles used was 105. The continuous dotted curve represents Eq. (61) plotted in the full range 0 < γ < 1. It holds, for example, for the indicator function, which is an integrable observable in this range, and hence, this observable does not exhibit the discontinuous behavior unlike the energy observable.

FIG. 9.

The EB parameter vs γ for the kinetic energy observable v2. When EB = 0, we have an ergodic phase found when γ > 1. This ergodicity breaking parameter is a measure for the fluctuations of the time averages, and for 1/3 < γ < 1, it is given by Eq. (61). That formula is valid whenever the observable is integrable with respect to the infinite measure. For 0 < γ < 1/3, the kinetic energy observable is non-integrable, and the system enters a different phase, described in Eq. (69). The dots are obtained from finite time simulations of a model with a parent velocity distribution f(v), which is uniform, and the number of particles used was 105. The continuous dotted curve represents Eq. (61) plotted in the full range 0 < γ < 1. It holds, for example, for the indicator function, which is an integrable observable in this range, and hence, this observable does not exhibit the discontinuous behavior unlike the energy observable.

Close modal

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

(70)

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̄=0tEk(t)dt/t=v2t/t=v2. Note that we consider here the limit where t is made long, and only then, γ → 0. It is easy to find ⟨Ek⟩ = ⟨v2⟩ and (Ek)2=v4 in the limit γ → 0 using Eqs. (19) and (20). In particular, exp(−x1/γ/c) = 1 if 0 < x < 1; otherwise, it is zero since γ → 0. Hence, using Eq. (19), v2(t)=01x0xdzdx(c/t)2γ=(1/3)(c/t)2γ and, similarly, ⟨v4(t)⟩ = (1/5)(c/t)4γ, and indeed, limγ→0EB = [(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.

We will now obtain the PDF of the random variable ϒ=S(t)/S(t)=Ēk(t)/Ēk(t), 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̂(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̂(u,p) attain scaling forms, denoted as P(S,t)t1α/2fα(S/t12/α) and P(u, p) ∼ (1/p)gα(u/p1−2/α). Here, the limit under study is t and S, the ratio S/t12/α remaining fixed and similarly in Laplace space. Note that we showed in Eq. (C8) that S(t)t12/α; 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,

(71)

Substituting x=St12/α and setting p = 1 turn Eq. (71) into the integral equation

(72)

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

(73)

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̂(u,p) [Eq. (53)] the following exact form:

(74)

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),

(75)

where

is the Mellin transform of gα(y) and Kα̃(s) and f̃α(s) is defined analogously. Solving Eq. (75) for f̃α(s) and applying the inverse Mellin transformation give

(76)

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

(77)

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

(78)

with the Mellin transform g̃4(s)=Γ(s)Γ(1s). Since the Mellin transform Kα̃(s) of the integral kernel is also known in full generality, K̃α(s)=Γ(s)Γ(1(12/α)s), we get the quotient in Eq. (76), and in addition, we can invert from Mellin space to obtain f4(x)=1πex24. The scaling factor C4=xf4 follows from the general relation between the nth derivative gα(n)(y=0) and the moments xnfα of fα(x),

(79)

which follows directly from Eq. (72). For α = 4, we get C4=xf4=2π so that for P4(ϒ), according to Eq. (77), we obtain a half Gaussian distribution

(80)

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

(81)

and eventually,

(82)

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, ϒ = v2/⟨v2⟩, 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.

FIG. 10.

Histogram of the normalized time averaged energy Pα(ϒ), obtained from numerical simulations with α = 50, is compared with the limiting PDF [Eq. (82)]. In the limit, the PDF is sharply cutoff at ϒ = 3, an effect which is smeared out with the finite time, finite α simulations. Note the low energy ϒ → 0 blow up of P(ϒ), a feature not found for α = 1.25 and α = 2 in Figs. 5 and 6 when the Mittag–Leffler law holds.

FIG. 10.

Histogram of the normalized time averaged energy Pα(ϒ), obtained from numerical simulations with α = 50, is compared with the limiting PDF [Eq. (82)]. In the limit, the PDF is sharply cutoff at ϒ = 3, an effect which is smeared out with the finite time, finite α simulations. Note the low energy ϒ → 0 blow up of P(ϒ), a feature not found for α = 1.25 and α = 2 in Figs. 5 and 6 when the Mittag–Leffler law holds.

Close modal

For also the experimentally relevant case α = 6, we can still get from Eq. (74) by the residue calculus a fully analytic expression for g6(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(y), which deviates only little from the true function g6(y), but for which the Mellin transform is known analytically. The general idea is to match in g6(y) the true small y-behavior and the large y-asymptotics of g6(y), which we know analytically from an analysis of Eq. (74). The simple form

(83)

shares with the exact function g6(y) the identical first and second derivative at y = 0 and the exponent of the asymptotic behavior for y. The relative deviation of g6(y) from g6(y) is strictly less than 4.2%. This function g6(y) can be Mellin transformed65 and leads via Eq. (76) and subsequent rescaling to an analytical expression for P6(ϒ), which can be expressed in terms of a Fox-H function66 as

(84)

As mentioned, after rewriting this expression in terms of a Meijer G-function (Ref. 65, p. 629), we can plot the result with Mathematica; see Fig. 11.

FIG. 11.

Simulations and theory Eq. (84) for the distribution of the time averaged energy Pα(ϒ) for α = 6. Here, the observable is non-integrable with respect to the infinite density, and P6(ϒ) diverges at ϒ → 0, indicating very long sticking times close to zero speed at least for some of the atoms that have very low energy compared to the mean.

FIG. 11.

Simulations and theory Eq. (84) for the distribution of the time averaged energy Pα(ϒ) for α = 6. Here, the observable is non-integrable with respect to the infinite density, and P6(ϒ) diverges at ϒ → 0, indicating very long sticking times close to zero speed at least for some of the atoms that have very low energy compared to the mean.

Close modal

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

(85)

with the decay exponent

(86)

and decay amplitude

(87)

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

(88)

with the exponent

(89)

and amplitude

(90)

The following table gives an overview for the resulting behavior of fα(x0)t(α)xs0(α) 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α(0),

(91)

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 P3+ɛ(x) as ɛ → 0 in contrast to the vanishing amplitude of f3+ɛ(x), which is a consequence of the diverging scaling factor C3+ɛɛ−1.

αfα(x → 0)Pα(x → 0)
α 12x12 123x12 
50 24sin(π48)πΓ(2548)x2348 24sin(π48)πΓ(2548)[(1+2cos(π25))Γ(4925)]2548x2348 
2πΓ(34)x14 1214πΓ(34)Γ(53)34x14 
334πΓ(56)x16 33[12(35)]5124πΓ(56)Γ(85)56x16 
1πx0 2πx0 
72 338πΓ(76)x16 338πΓ(76)cos(3π14)sin(π7)Γ(107)76x16 
103 23πΓ(54)x14 (1+5)5212πΓ(54)Γ(75)54x14 
α = 3 + ɛ εx12ε (23πΓ(75))32ε12x12ε 
αfα(x → 0)Pα(x → 0)
α 12x12 123x12 
50 24sin(π48)πΓ(2548)x2348 24sin(π48)πΓ(2548)[(1+2cos(π25))Γ(4925)]2548x2348 
2πΓ(34)x14 1214πΓ(34)Γ(53)34x14 
334πΓ(56)x16 33[12(35)]5124πΓ(56)Γ(85)56x16 
1πx0 2πx0 
72 338πΓ(76)x16 338πΓ(76)cos(3π14)sin(π7)Γ(107)76x16 
103 23πΓ(54)x14 (1+5)5212πΓ(54)Γ(75)54x14 
α = 3 + ɛ εx12ε (23πΓ(75))32ε12x12ε 

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.

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/Tg]/Tg, where E > 0 are the trap depths and Tg 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=0ρ(E)exp(E/T)dE (here, kB = 1), and hence, it diverges when Tg > 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.

  1. It is a challenge to see if quantum Monte Carlo simulations7 can be used to numerically investigate the non-normalizable state and the time and ensemble averages.

  2. Experimentally finding the infinite density is another obvious challenge. Single atom experiments yield direct insights into the trajectories and the time averages.68–71 

  3. 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 considered55 other models, including the case where the process is not renewed after each jolt. The main results are left unchanged.

  4. 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?

  5. 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 ta, and only then perform a measurement, i.e., time average in the interval (ta, ta + t). In this case, we expect that statistical properties of the system will depend on the aging time ta. 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 

  6. 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.

  7. 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.

  8. 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(va < v(t) < vb) is integrable with respect to the infinite density when 0 < va and non-integrable if va = 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(−kx2/2kBT); 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.

  9. 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 vd−1 for v < vmax 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)vα+d1, 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.

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̄(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.

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.

The authors have no conflicts to disclose.

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

1. Basics of renewal theory

We start with a brief recapitulation of renewal theory.44,45,73–75 Let ϕ1(τ̃) 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 QN(t)dt be the probability that the Nth renewal event is taking place in the interval (t, t + dt). Then, from the renewal property of the process, we have

(A1)

with the initial condition Q0(t) = δ(t). The probability of finding N renewals in (0, t) is

(A2)

Here, W(t)=10tϕ1(τ̃) is the probability of not making a transition up to time t. Equation (A2), thus, describes a situation where the Nth renewal takes place at time tτ̃, and in the remaining time τ̃, no jump was made. We now consider the Laplace transform of PN(t) denoted as P̂N(p), and using Eqs. (A1) and (A2) and the convolution theorem, we find

(A3)

It follows that the mean of N is

(A4)

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

(A5)

where in our case, τ0=c[Γ(1+1/α)/vmax]α. Then, in the limit p → 0, one can show that for α > 1,

(A6)

where the leading term is the normalization. Inserting Eq. (A6) into Eq. (A4) and then performing a straightforward inverse Laplace transform, one finds N(t)[αsin(π/α)/π](t/τ0)1/α.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 ϕ̂1(p)1bγpγ and, then,

(A7)

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

(A8)

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

(A9)

and using convolution, Φ̂2(0,p)=[1ϕ̂2(0,p)]/p. Hence, from Eq. (55), we have

(A10)

When u = 0, the Laplace transform of the joint PDF ϕ̂2(u,p) reduced to the Laplace transform of the marginal PDF of waiting times, namely, ϕ̂2(0,p)=ϕ̂1(p). We now use Eq. (A4) and ϕ̂1(p)1 to leading order in p, finding Ŝ(p)sN̂(p), which gives Eq. (56).

We consider the fluctuations of the time average of the energy for 1 < α < 3. Since Ek̄=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,

(B1)

Integrating over s and τ̃ gives

(B2)

Similarly, we find

(B3)

For u = 0, we have Φ̂2(0,p)=[1ϕ̂2(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γ, 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

(B4)

As communicated, we consider the case 1/3 < γ < 1, and in this regime, p2γ−1u = p3γ−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

(B5)

and 0(1+zγ)1dz=γπcsc(γπ) with γ < 1.

We now consider the denominator of Eq. (53) and note that limu,p0ϕ̂2(u,p)=1 from normalization. Rewriting, we find

(B6)

In the limit G1(cγpγ/vmax)0dz(1+zγ)1 and G2c(vmax)2αu/(3α) or using Eq. (44), G2us⟩. We see that G1pγ and G2u, so these two terms are of the same order.

Inserting Eqs. (B5) and (B6) into Eq. (53), we find

(B7)

with (T1)γ=γπcγ/(vmaxsinγπ). Now, the inverse Laplace transform uS gives

(B8)

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/dpt, we obtain the PDF of the action for 1/3 < γ < 1, sometimes called the Mittag–Leffler law,

(B9)

with t̃=t/T1 and S̃=S/s. 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.

1. Formula for Ēk

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 Ŝ(p), 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 Ēk(t)S(t)/t.

Equation (55) is split into two terms,

(C1)

with

(C2)

When 1/3 < γ, the derivative uϕ̂2(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 Φ̂2(u,p) is associated with the probability of not jolting, namely, the term Ĥ2(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),

(C3)

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

(C4)

The upper limit is then taken to be infinity, and we find in the small p limit ϕ̂2(0,p)1(pc)γπγ/vmaxsinπγ. 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

(C5)

all valid in the regime of interest.

Using the same tricks, we find

(C6)

and

(C7)

Inserting Eqs. (C6) and (C8) into Eq. (C2), we find the small p limit of Ĥ1(p) and Ĥ2(p). It is now easy to transform from p to time t and find

(C8)

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 time53,54 [the H2(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 Ij(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

(C9)

and using Eq. (C4),

(C10)

Note that I1(p) is independent of vmax in this limit, and so are the remaining terms I2(p), I3(p) and I4(p), which are given by

(C11)
(C12)

and similarly,

(C13)

Inverting to the time domain, we see that ⟨S2(t)⟩∝ t2−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 ⟨S2(t)⟩ − ⟨S(t)⟩2, and this after normalization yields the max EB parameter [Eq. (69)].

1.
S.
Chu
, “
The manipulation of neutral particles
,”
Rev. Mod. Phys.
70
,
685
(
1998
).
2.
C. N.
Cohen-Tannoudji
, “
Manipulating atoms with photons
,”
Rev. Mod. Phys.
70
,
707
(
1998
).
3.
W. D.
Phillips
, “
Laser cooling and trapping of neutral atoms
,”
Rev. Mod. Phys.
70
,
721
(
1998
).
4.
E. S.
Shuman
,
J. F.
Barry
, and
D.
DeMille
, “
Laser cooling of a diatomic molecule
,”
Nature
467
,
820
(
2010
).
5.
F.
Bardou
,
J. P.
Bouchaud
,
O.
Emile
,
A.
Aspect
, and
C.
Cohen-Tannoudji
, “
Sub-recoil laser cooling and Lévy flights
,”
Phys. Rev. Lett.
72
,
203
(
1994
).
6.
S.
Marksteiner
,
K.
Ellinger
, and
P.
Zoller
, “
Anomalous diffusion and Lévy walks in optical lattices
,”
Phys. Rev. A
53
,
3409
(
1996
).
7.
F.
Bardou
,
J. P.
Bouchaud
,
A.
Aspect
, and
C.
Cohen-Tannoudji
,
Lévy Statistics and Laser Cooling: How Rare Events Bring Atoms to Rest
(
Cambridge University Press
,
2002
).
8.
E.
Lutz
and
F.
Renzoni
, “
Beyond Boltzmann–Gibbs statistical mechanics in optical lattices
,”
Nat. Phys.
9
,
615
(
2013
).
9.
G.
Afek
,
N.
Davidson
,
D. A.
Kessler
, and
E.
Barkai
, “
Anomalous statistics of laser-cooled atoms in dissipative optical lattices
,” arXiv:2107.09526 [cond-mat.stat-mech].
10.
D. A.
Darling
and
M.
Kac
, “
On occupation times for Markoff process
,”
Trans. Am. Math. Soc.
84
,
444
(
1957
).
11.
J.
Aaronson
,
An Introduction to Infinite Ergodic Theory
(
AMS
,
1997
).
12.
R.
Zweimüller
,
Surrey Notes on Infinite Ergodic Theory
, Lecture Notes (
University of Surrey
,
2009
).
13.
N.
Korabel
and
E.
Barkai
, “
Pesin-type identity for intermittent dynamics with a zero Lyapunov exponent
,”
Phys. Rev. Lett.
102
,
050601
(
2009
).
14.
D. A.
Kessler
and
E.
Barkai
, “
Infinite covariant density for diffusion in logarithmic potentials and optical lattices
,”
Phys. Rev. Lett.
105
,
120602
(
2010
).
15.
T.
Akimoto
and
T.
Miyaguchi
, “
Role of infinite invariant measure in deterministic sub diffusion
,”
Phys. Rev. E
82
,
030102
(
2010
).
16.
T.
Akimoto
, “
Distributional response to biases in deterministic superdiffusion
,”
Phys. Rev. Lett.
108
,
164101
(
2012
).
17.
P.
Meyer
and
H.
Kantz
, “
Infinite invariant densities due to intermittency in a nonlinear oscillator
,”
Phys. Rev. E
96
,
022217
(
2017
).
18.
A.
Vezzani
,
E.
Barkai
, and
R.
Burioni
, “
Single-big-jump principle in physical modeling
,”
Phys. Rev. E
100
,
012108
(
2019
).
19.
E.
Aghion
,
D. A.
Kessler
, and
E.
Barkai
, “
From non-normalizable Boltzmann–Gibbs statistics to infinite-ergodic theory
,”
Phys. Rev. Lett.
122
,
010601
(
2019
);
[PubMed]
E.
Aghion
,
D. A.
Kessler
, and
E.
Barkai
Infinite ergodic theory meets Boltzmann statistics
,”
Chaos, Solitons Fractals
138
,
109890
(
2020
).
20.
Y.
Sato
and
R.
Klages
, “
Anomalous diffusion in random dynamical systems
,”
Phys. Rev. Lett.
122
,
174101
(
2019
).
21.
T.
Akimoto
,
E.
Barkai
, and
G.
Radons
, “
Infinite invariant density in a semi-Markov process with continuous state variables
,”
Phys. Rev. E
101
,
052112
(
2020
).
22.
M.
Radice
,
M.
Onofri
,
R.
Artuso
, and
G.
Pozzoli
, “
Statistics of occupation times and connection to local properties of non-homogeneous random walks
,”
Phys. Rev. E
101
,
042103
(
2020
).
23.
I.
Fouxon
and
P.
Ditlevsen
, “
Refined central limit theorem and infinite density tail of the Lorentz gas from Lévy walk
,”
J. Phys. A
53
,
415004
(
2020
).
24.
J. P.
Bouchaud
, “
Weak ergodicity breaking and aging in disordered systems
,”
J. Phys. I: EDP Sci.
2
(
9
),
1705
1713
(
1992
).
25.
R.
Metzler
,
J.-H.
Jeon
,
A. G.
Cherstvy
, and
E.
Barkai
, “
Anomalous diffusion models and their properties: Non-stationarity, non-ergodicity and ageing at the centenary of single particle tracking
,”
Phys. Chem. Chem. Phys.
16
(
44
),
24128
24164
(
2014
).
26.
E.
Bertin
and
F.
Bardou
, “
From laser cooling to aging: A unified Lévy flight description
,”
Am. J. Phys.
76
,
630
(
2008
).
27.
J.
Reichel
et al., “
Raman cooling of cesium below 3 nK: New approach inspired by Lévy flight statistics
,”
Phys. Rev. Lett.
75
,
4575
(
1995
).
28.
J.
Reichel
, “
Refroidissement Raman et vols de Lévy: Atomes de césium au nanoKelvin
,” Physique Atomique [physics.atom-ph], Université Pierre et Marie Curie—Paris VI, 1996, Francais, Tel.:00004691.
29.
W. E.
Moerner
and
M.
Orrit
, “
Illuminating single molecules in condensed matter
,”
Science
283
,
1670
(
1999
).
30.
E.
Barkai
,
Y.
Garini
, and
R.
Metzler
, “
Strange kinetics of single molecules in the cell
,”
Phys. Today
65
(
8
),
29
(
2012
).
31.
P.
Frantsuzov
,
M.
Kuno
,
B.
Jankó
, and
R. A.
Marcus
, “
Universal emission intermittency in quantum dots, nanorods and nanowires
,”
Nat. Phys.
4
,
519
(
2008
).
32.
F. D.
Stefani
,
J. P.
Hoogenboom
, and
E.
Barkai
, “
Beyond quantum jumps: Blinking nano-scale light emitters
,”
Phys. Today
62
(
2
),
34
(
2009
).
33.
E.
Barkai
,
G.
Radons
, and
T.
Akimoto
, “
Transitions in the ergodicity of subrecoil-laser-cooled gases
,”
Phys. Rev. Lett.
127
,
140605
(
2021
).
34.
G. D.
Birkhoff
, “
Proof of the ergodic theorem
,”
Proc. Natl. Acad. Sci. U. S. A.
17
,
656
(
1931
).
35.
D.
Dijkstra
, “
A continued fraction expansion for a generalisation of Dawson’s integral
,”
Math. Comput.
31
,
503
(
1977
).
36.
I. S.
Gradshteyn
and
I. M.
Ryzhik
,
Table of Integrals, Series, and Products
, 7th ed. (
Elsevier Academic Press
,
Amsterdam
,
2007
).
37.
C.
Cohen-Tannoudji
and
D.
Guéry-Odelin
, in
Advances in Atomic Physics an Overview
(
World Scientific, Singapore
,
2011
), pp.
291
315
.
38.
B.
Saubaméa
,
M.
Leluc
, and
C.
Cohen-Tannoudji
, “
Experimental investigation of non-ergodic effects in sub recoil laser cooling
,”
Phys. Rev. Lett.
83
,
3796
(
1999
).
39.
R.
Metzler
and
J.
Klafter
, “
Random walk’s guide to anomalous diffusion: A fractional dynamics
,”
Phys. Rep.
339
,
1
(
2000
).
40.
R.
Kutner
and
J.
Masoliver
, “
The continuous time random walk, still trendy: Fifty-year history, state of art and outlook
,”
Eur. Phys. J. B
90
,
50
(
2017
).
41.
V.
Zaburdaev
,
S.
Denisov
, and
J.
Klafter
, “
Lévy walks
,”
Rev. Mod. Phys.
87
,
483
(
2017
).
42.
T.
Akimoto
,
S.
Shinkai
, and
Y.
Aizawa
, “
Distributional behavior of time averages of non-L1 observables in one-dimensional intermittent maps with infinite invariant measures
,”
J. Stat. Phys.
158
,
476
(
2015
).
43.
T.
Albers
and
G.
Radons
, “
Exact results for the nonergodicty of d-dimensional generalized Lévy walks
,”
Phys. Rev. Lett.
120
,
104501
(
2018
).
44.
C.
Godrèche
and
J. M.
Luck
, “
Statistics of the occupation time of renewal processes
,”
J. Stat. Phys.
104
,
489
(
2001
).
45.
W.
Wang
,
J. H. P.
Schulz
,
W.
Deng
, and
E.
Barkai
, “
Renewal theory with fat tailed distributed sojourn times: Typical versus rare
,”
Phys. Rev. E
98
,
042139
(
2018
).
46.
J.
Klafter
,
A.
Blumen
, and
M. F.
Shlesinger
, “
Stochastic pathway to anomalous diffusion
,”
Phys. Rev. A
35
,
3081
(
1987
).
47.
T.
Akimoto
and
T.
Miyaguchi
, “
Distributional ergodicity in stored-energy-driven Lévy flights
,”
Phys. Rev. E.
87
,
062134
(
2013
).
48.
T.
Akimoto
and
T.
Miyaguchi
, “
Phase diagram in stored-energy-driven Lévy flight
,”
J. Stat. Phys.
157
,
515
(
2014
).
49.
E.
Aghion
,
D. A.
Kessler
, and
E.
Barkai
, “
Asymptotic densities from the modified Montroll–Weiss equation for coupled CTRWs
,”
Eur. Phys. J. B
91
,
17
(
2018
).
50.
E. W.
Montroll
and
G. H.
Weiss
, “
Random walks on lattices. II
,”
J. Math. Phys.
6
,
167
(
1965
).
51.
K. A.
Penson
and
K.
Górska
, “
Exact and explicit probability densities for one-sided Lévy stable distributions
,”
Phys. Rev. Lett.
105
,
210604
(
2010
).
52.
Y.
He
,
S.
Burov
,
R.
Metzler
, and
E.
Barkai
, “
Random time-scale invariant diffusion and transport coefficients
,”
Phys. Rev. Lett.
101
,
058101
(
2008
).
53.
C.
Godrèche
,
S. N.
Majumdar
, and
G.
Schehr
,
J. Stat. Mech.
2015
,
P03014
.
54.
M.
Höll
,
W.
Wang
, and
E.
Barkai
, “
Extreme value statistics for constrained physical models
,”
Phys. Rev. E
102
,
042141
(
2020
).
55.
T.
Akimoto
,
E.
Barkai
, and
G.
Radons
, “
Infinite ergodic theory for three heterogeneous stochastic models with application to subrecoil laser cooling
,” (unpublished).
56.
Y.
Castin
,
J.
Dalibard
, and
C.
Cohen-Tannoudji
, “
The limits of Sisyphus cooling
,” in
Light Induced Kinetic Effects on Atoms, Ions and Molecules
, edited by
L.
Moi
, et al.
(
ETS Editrice
,
Pisa
,
1991
).
57.
P.
Douglas
,
S.
Bergamini
, and
F.
Renzoni
, “
Tunable tsallis distributions in dissipative optical lattices
,”
Phys. Rev. Lett.
96
,
110601
(
2006
).
58.
Y.
Sagi
,
M.
Brook
,
I.
Almog
, and
N.
Davidson
, “
Observation of anomalous diffusion and fractional self-similarity in one dimension
,”
Phys. Rev. Lett.
108
,
093002
(
2012
).
59.
D. A.
Kessler
and
E.
Barkai
, “
Theory of fractional-Lévy kinetics for cold atoms diffusing in optical lattices
,”
Phys. Rev. Lett.
108
,
230602
(
2012
).
60.
A.
Dechant
and
E.
Lutz
, “
Anomalous spatial diffusion and multifractality in optical lattices
,”
Phys. Rev. Lett.
108
,
230601
(
2012
).
61.
E.
Barkai
,
E.
Aghion
, and
D.
Kessler
, “
From the area under the Bessel excursion to anomalous diffusion of cold atoms
,”
Phys. Rev. X
4
,
021036
(
2014
).
62.
P. C.
Holz
,
A.
Dechant
, and
E.
Lutz
, “
Infinite density for cold atoms in shallow optical lattices
,”
Europhys. Lett.
109
,
23001
(
2015
).
63.
T.
Akimoto
and
E.
Barkai
, “
Aging generates regular motions in weakly chaotic systems
,”
Phys. Rev. E.
87
,
032915
(
2013
).
64.
A. D.
Polyanin
and
A. V.
Manzhirov
,
Handbook of Integral Equations
, 2nd ed. (
Chapman & Hall; CRC
,
Boca Raton
,
2008
).
65.
A. P.
Prudnikov
,
Yu. A.
Brychkov
, and
O. I.
Marichev
,
Integrals and Series
, More Special Functions Vol. 3 (
Gordon & Breach
,
New York; London
,
1989
).
66.
A. M.
Mathai
,
R. K.
Saxena
, and
H. J.
Haubold
,
The H-Function Theory and Applications
(
Springer
,
New York
,
2010
).
67.
G.
Bel
and
E.
Barkai
, “
Weak ergodicity breaking in the continuous-time random walk
,”
Phys. Rev. Lett.
94
,
240602
(
2005
).
68.
H.
Katori
,
S.
Schlipf
, and
H.
Walther
, “
Anomalous dynamics of a single ion in an optical lattice
,”
Phys. Rev. Lett.
79
,
2221
(
1997
).
69.
I.
Stroescu
,
D. B.
Hume
, and
M. K.
Oberthaler
, “
Dissipative double-well potential for cold atoms: Kramers rate and stochastic resonance
,”
Phys. Rev. Lett.
117
,
243005
(
2016
).
70.
M.
Hohmann
 et al., “
Single-atom thermometer for ultracold gases
,”
Phys. Rev. A
93
,
043607
(
2016
).
71.
Q.
Bouton
 et al., “
Single-atom quantum probes for ultracold gases boosted by nonequilibrium spin dynamics
,”
Phys. Rev. X
10
,
011018
(
2020
).
72.
N.
Davidson
,
H. J.
Lee
,
M.
Kasevich
, and
S.
Chu
, “
Raman cooling of atoms in two and three dimensions
,”
Phys. Rev. Lett.
72
,
3158
(
1994
).
73.
D. R.
Cox
,
Renewal Theory
(
Methuen and Co., Ltd.
,
London
,
1962
).
74.
J. E.
Santos
,
T.
Franosch
,
A.
Parmeggiani
, and
E.
Frey
, “
Renewal processes and fluctuation analysis of molecular motor stepping
,”
Phys. Biol.
2
,
207
222
(
2005
).
75.
D.
Daley
and
D.
Vere-Jones
,
An Introduction to the Theory of Point Processes
, Elementary Theory and Methods Vol. I (
Springer
,
2003
).