Infinite ergodic theory meets Boltzmann statistics

Chaos, Solitons and Fractals 138 (2020) 109890 Contents lists available at ScienceDirect Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos Frontiers Infinite ergodic theory meets Boltzmann statistics Erez Aghion a,b,c,∗, David A. Kessler a, Eli Barkai a,b a Department of Physics, Institute of Nanotechnology and Advanced Materials, Bar-Ilan University, Ramat-Gan 52900, Israel Institute of Nanotechnology and Advanced Materials, Bar-Ilan University, Ramat-Gan 52900, Israel c Max-Planck Institute for the Physics of Complex Systems, Dresden D-01187, Germany b a r t i c l e i n f o Article history: Received 4 February 2020 Revised 5 May 2020 Accepted 10 May 2020 a b s t r a c t We investigate the overdamped stochastic dynamics of a particle in an asymptotically flat external potential field, in contact with a thermal bath. For an infinite system size, the particles may escape the force field and diffuse freely at large length scales. The partition function diverges and hence the standard canonical ensemble fails. This is replaced with tools stemming from infinite ergodic theory. BoltzmannGibbs statistics, even though not normalized, still describes integrable observables, like energy and occupation times. The Boltzmann infinite density is derived heuristically using an entropy maximization principle, as well as via a first-principles calculation using an eigenfunction expansion in the continuum of low-energy states. A generalized virial theorem is derived, showing how the virial coefficient describes the delay in the diffusive spreading of the particles, found at large distances. When the process is non-recurrent, e.g. diffusion in three dimensions with a Coulomb-like potential, we use weighted time averages to restore basic canonical relations between time and ensemble averages. © 2020 Elsevier Ltd. All rights reserved. 1. Introduction The overdamped stochastic dynamics of a particle in an external potential field V(x), in one dimension, in contact with a thermal bath, is given by the Langevin equation V ′ (x ) √ dx =− + 2Dη (t ). dt γ (1) Here −V ′ (x ) is the deterministic force applied on the particle due to the potential field, and γ > 0 is the friction constant. η(t) is the bath noise, which is white, Gaussian, has zero mean and η (t )η (t ′ ) = δ (t − t ′ ) (where δ ( · ) is the Dirac δ -function). The Einstein relation D = kB T /γ guarantees that the system, in the case of a binding potential V(x), will relax to thermal equilibrium. In this case, the steady-state equilibrium density is [1] Peq (x ) = exp [−V (x )/kB T ] . Z (2) This final equilibrium state transcends a particular type of dynamics, and the asymptotic shape of the density does not depend on transport coefficients, such as the diffusion constant D, of the particles in the medium. Here, Z= ∗  ∞ −∞ exp [−V (x )/kB T ]dx Corresponding author. E-mail address: [email protected] (E. Aghion). https://doi.org/10.1016/j.chaos.2020.109890 0960-0779/© 2020 Elsevier Ltd. All rights reserved. (3) is the normalizing partition function, and kB is the Boltzmann constant. A finite value of Z is, however, not always guaranteed. In particular, Z can diverge when V(x) generates a force F (x ) = −V ′ (x ) → 0 when x → ∞ and/or −∞. More specifically, in this manuscript we are interested in the case where V(x) itself drops to zero at large distance, at least as fast as 1/x. We initiated a study of this case in a previous work [2], finding that at long times, and finite x, the time-dependent density assumes the shape of the BoltzmannGibbs factor, multiplied by a factor which decays as power-law in time. In the limit t → ∞, the Boltzmann-Gibbs factor becomes an infinite-invariant density [2]. In the potential-free region, the density is simply the free-diffusion kernel. The appearance of the Boltzmann-Gibbs density can be understood as resulting from the fact that the particle returns infinitely many times to the potential region and so a kind of conditioned equilibrium is established there. In addition, we then showed in [2] how in this thermal setting, we recover the Aaronson-Darling-Kac theorem [3,4], which describes the ergodic properties of a certain class of observables, integrable with respect to the infinite density. In the current work, we extend our previous study in several directions. Most notably, we re-derive our previous (one-dimensional) results using an eigenfunction expansion of the relevant time-dependent FokkerPlanck equation and thereby not only succeed in recovering the infinite-invariant density, but the leading-order corrections as well. We also derive a virial theorem for our system, and extend some of our results to d ≥ 1 dimensions. 2 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 1. Various classes of potential landscapes: a. An example of a binding potential; V(x) ∼ x2 . Here, the standard Boltzmann-Gibbs equilibrium state is achieved in the longtime limit, and standard ergodic theory applies. b. A logarithmic potential, V(x) ∼ log (|x|) for |x| ≫ 1; here, the Boltzmann-Gibbs distribution is normalizable if the temperature is sufficiently small, otherwise it leads to non-normalizable Boltzmann-Gibbs statistics [2,5,6]. c.+d. Two-sided and one-sided asymptotically flat potentials always lead to a non-normalizable Boltzmann-Gibbs state, Eq. (24), and infinite-ergodic theory (see also [2,7]). e. A periodic potential, with a structure that stretches for x ∈ (−∞, ∞ ) (also here one finds a non-normalised state, see e.g., [8]). f. An unstable potential (see e.g., [9,10]). All the potentials c-e, share the common property that in one dimension the process x(t) in Eq. (1) is recurrent, however the mean return-time from large x into any finite region around the origin is infinity. Thus, these potentials are weakly binding. This a sufficient, albeit not necessary condition for the emergence of an infinite-invariant density. There are many other situations where Z diverges as well. Some examples are presented in Fig. 1, to be compared with the binding potential Fig. 1a which leads to finite Z. For logarithmic potentials, such as the example in Fig. 1b, the partition function diverges when the depth of the well is sufficiently shallow. This happens when V(x) ∼ V0 log (|x|) at large |x|, and V0 < kB T at a certain given temperature, and the infinite invariant density of this class of potentials was studied in [2,5,6]). In addition, one may consider other non-binding fields, for example periodic structures, Fig. 1e, random potentials, or unstable fields, Fig. 1f. All the examples ce share two properties: first, the partition function diverges, and second, the Langevin dynamics is recurrent in one dimension. In turn, in one dimension this implies that the mean return time in these examples is diverging. We will refer to the class of potentials which fulfill these conditions as weakly binding. Logarithmic potentials are a marginal case, which behaves as binding or weakly binding given the system parameters, and the unstable potential in Fig. 1f belongs to neither group. As will become clear below, while binding potentials lead to standard ergodic theory, we anticipate that infinite ergodic theory will serve as a useful tool for Langevin dynamics in one dimension for all the weakly binding potentials (though as mentioned, in this work we study in detail only asymptotically flat fields). Note that the extension of the theory to higher dimensions is not only a technical issue, since in that case the random walk is no longer recurrent. The basic tools to deal with this case need some modifications, as we demonstrate for isotropic potentials whose amplitude drops to zero at large radial distances r, at least as fast as 1/r, below. The structure of the manuscript is as follows: We first discuss in Section 2 some preliminary matters, presenting a brief recap of equilibrium statistical mechanics for binding potentials, where the partition function is normalizable, together with a description of the particular examples of asymptotically flat potentials we use for our simulations. In Section 3 we define the nonnormalizable Boltzmann state. In Section 4, we discuss the entropy maximization principle. In Section 5, we provide the derivation of the leading-order time-dependent shape of the particle density and obtain higher-order correction terms. In Section 6, we discuss time and ensemble averages, and in Section 7 we discuss the fluctuations of the time average and infinite ergodic theory. In Section 8, we study the virial theorem. In Section 9, we show that the non-normalizable Boltzmann state exists in any dimension d ≥ 1, and extend our analysis of the ratio between time and ensemble averages of integrable observables, to any dimension. A summary of our main results is found in Section 10. The discussion is found in Section 11. 2. Preliminaries 2.1. A recap of statistical mechanics Before treating weakly binding potentials, we first recall the standard treatment of the case where V(x) is increasing with distance in such a way that the partition function is finite, e.g., V (x ) = x2 (a binding potential, Fig. 1.a). According to the basic laws of statistical physics, the system is ergodic (we assume that V(x) does not divide the system into compartments). Let O[x(t )] be a physical observable. Then, in the long-time limit, O[x(t )] → O (x ). (4) O[x(t )] = Here, the overline denotes a time-average  1 t ′ )]dt ′ , and the brackets an ensemble-average. The O[ x ( t t 0 fact that the time-average, which is what is measured in many experiments, converges to the corresponding ensemble-average (in the long-time limit), is very useful for the theoretician, who usually considers the latter; O (x ) = ∞ 0 O (x ) exp[−V (x )/kB T ]dx . Z (5) It should be noted that some observables, such as O (x ) = exp[V (x )/T ], are not integrable with respect to the Boltzmann distribution, however these are mostly not the main focus of physicists. Statistical mechanics is related to thermodynamics in many textbooks. In particular, the Helmholtz free energy is (6) F = V  − T S, where we omitted the thermal kinetic energy term without any loss of generality. The entropy is S(t ) = −kB  ∞ 0 Pt (x ) ln[Pt (x )]dx, (7) where Pt (x)dx is the probability of finding the particle at time t in the interval (x, x + dx ). In equilibrium, we take the long-time limit, and for generic initial conditions we have lim Pt (x ) = Peq (x ), t→∞ (8) so, in this limit, F = −kB T ln Z, (9) and Peq (x ) = e−V (x )/kb T . e−F/kb T (10) The Boltzmann factor appearing in the numerator is the essence of the canonical ensemble, while in the denominator we have the relation to thermodynamics. The goal of this manuscript is to consider the case where S is increasing with time, as opposed to saturating to a limit, but still all this structure remains intact when the appropriate modifications are made, namely we must use the tool of non-normalizable Boltzmann-Gibbs statistics [2,5,7]. This idea, which is discussed at length below, harnesses the tools of infinite-ergodic theory, which has been well established as a mathematical theory for several decades, and was discovered in recent years also in physical systems, e.g., [3,4,11–22]. Note that infinite, namely non-normalizable, densities serve two main goals; the first is for computation of large deviations and rare-event statistics of fat-tailed stochastic systems, and the second is in the context of E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 3 ergodic theory. In that regard, one should distinguish between infinite invariant and infinite covariant densities, the latter are not discussed here (see e.g., Refs. [23–30]). 2.2. Asymptotically flat potentials As mentioned above, in this work we treat the class of potentials where V(x) → 0 at x → ∞, and the drop rate is at least as fast as 1/x. In this case V(x) belongs to the larger class of potentials which are weakly binding. Our leading-order results below apply to any such asymptotically flat potential; however, some of our calculations of higher-order correction terms below are obtained only for potentials that fall-off faster than 1/x2 . Furthermore, we distinguish between two situations: the one-sided case, where V (x ) = ∞ when x ≤ 0, and the double-sided system, where, to simplify explanations, we assume that the system is symmetric, and so V(x) → 0 also when x → −∞ (but our results can be trivially extended also to the non-symmetric case). The first situation can be realized experimentally, for example, when the particles are diffusing in three dimensions in a heat bath, above a flat surface, and their interaction with that surface is given by V(x), where x is their height (see e.g., [31–36]). Here, the potential is infinite when x is zero, since the surface is impenetrable. The second case, corresponds e.g., to a scenario where the particles diffuse in a liquid while being loosely held by optical tweezers. The intensity of such an optical trap drops with the distance from the potential well, often like an inverted Gaussian [37–39]. In our simulations, we mostly used the following two examples: the one-sided Lennard-Jones type potential (which depends only on the height coordinate, x, in the scenario of three dimensional diffusion above a hard-surface): V (x ) =  V0  12∞  6  a b x − x x≤0 x>0 , (11) where V0 , a and b are positive constants, and the symmetric potential (which can be realized using optical tweezers e.g., [38,39]): 2 V (x ) = (x/5 )4 − (x/5 )2 e−(x/5) . (12) In both cases, there is no thermal equilibrium in the usual sense. Famously, this “problem” with asymptotically flat potentials was pointed out by Fermi already in 1924 [40]. Physically, when x is large, the deterministic force is negligible and then the particles are diffusing in the bulk. Our discussion below is not limited to a specific form of the potential field, provided that it is eventually flat, but the key assumption is that the fluctuationdissipation relation applies, namely the Einstein relation is valid. This fact, as we show, allows for modified thermal concepts to emerge, even though Z = ∞. Note that unless specified explicitly, below we present our derivations mostly for the one-sided potentials, just for simplicity of writing. 2.3. Limitations of the standard treatment of the non-normalizability of the partition function, and the alternative The standard response to the “problem” of the divergence of the Boltzmann factor, for any type of potential, is to introduce a finite size to the system, L, and since then the limits of the integral in Eq. (3) will stretch only up to this limit, it is guaranteed to have a finite value. In this case, the system will always approach equilibrium in the long-time limit. On the bottom-left panel of Fig. 2, we see an example of a system enclosed between two walls, with a potential that, without the walls (namely, if the system had stretched from −∞ to ∞), would have been weakly binding. On the top-left panel of the figure, we see the corresponding Boltzmann distribution of the particles. But in this work, we Fig. 2. In an open system, stretching from −∞ to ∞, the potential V (x ) = [(x/5 )4 − 2(x/5 )2 ] exp[−(x/5 )2 ] is weakly binding, like the one seen in Fig. 1c. There are two ways by which one can obtain the Boltzmann distribution in a system with such a potential. Left bottom corner: we place the system within hard walls at ± 10, and since the system has a finite size, it has to eventually relax to equilibrium. Upper left corner: we see the simulated position distribution of the particles (blue circles), where the normalization is obtained by dividing the histogram of the positions by the number of particles and a bin-size. The theoretical red curve corresponds to Eq. (2). Right bottom corner: when we remove the hard walls, particles escape freely to infinity. This case is very different from the previous one as here the system size is infinite. However, here we track only the section of −10 < x < 10, and the particles which are found there at time t. Upper right corner: The simulated position distribution of the particles, normalized only with respect to the particles that are found within the “illuminated” region of x ∈ [−10, 10] (blue squares) at time t = 10 0 0 0, versus the same theoretical red curve as above. As predicted by our previous results in Ref. [2], the upper two pictures are practically indistinguishable, despite describing very different physical scenarios. Our goal in this paper is to investigate in depth the dynamics of the latter scenario. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) do not wish to impose this constraint. One reason is simply that many experimental settings do not have truly hard walls or confining potentials e.g., when using optical traps (see for example, [38,39]). The second reason is that even when the systems is confined by its boundaries, very often particle tracking experiments are not long enough for the particles to encounter them. The results of our previous work, Ref. [2], imply that we can take a very different physical approach, but still retain the exact same shape of Pt (x) at long times, obtained in the hard-wall scenario. In this second approach, we focus our range of observation to a single slice of space, and normalize the probability distribution only with respect to the number of particles that are present in this region at time t (namely, by discarding the particles which are found outside of this region). Note, that this situation is common in single-particle-tracking experiments, where the microscope in use has a finite field of view, hence this approach is essentially similar to “looking under the lamp”. The result of a simulation of this is shown in the right side of Fig. 2, to be compared with the wall scenario on the left. The measured concentration of the remaining particles is identical (up to statistical fluctuations) to that of the equilibrium state of the “walled” system, even though the distributions come from two different physical setups. As we shall see in the following, the dynamics by which the two final distributions are attained are extremely different. The “walled” system converges exponentially in time (with the time diverging as L → ∞), whereas the alternative converges as a power-law. 3. The non-normalized boltzmann-Gibbs state In this section, we review the derivation of the long-time limit of the distribution function to leading-order. We then analyze its thermodynamic implications. We start from the Fokker-Planck equation description of the diffusion process controlled by the Langevin Eq. (1), which specifies the dynamics of the concentra- 4 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 tion of particles, or equivalently the probability density function Pt (x), ∂2 ∂ V ′ (x ) ∂ Pt (x ) =D P ( x ). + 2 ∂t ∂ x kB T t ∂x (13) We treat this equation in the long-time limit. If we set the lefthand side to zero, namely we search for a time-independent solution, which we call I (x ), we have 0= ∂ V ′ (x ) ∂2 T I ( x ), + 2 ∂ x kB ∂x (14) and hence one appealing option reads I (x ) = Const exp [−V (x )/kB T ]. (15) However, this solution does not satisfy the boundary condition limx→∞ Pt (x ) = 0 (unless Const = 0) and is clearly nonnormalizable when V(x) is asymptotically flat at large distances. This is certainly not a possibility as the particles are neither created nor annihilated, so the normalization is conserved ∞ 0 Pt (x )dx = 1 for any t ≥ 0. In-fact, mathematically as we will show below, the non-normalized solution I (x ) is an infiniteinvariant density [4], as opposed to a probability density. Since I (x ) is non-normalized, we search for a more complete solution in the form of − V (x ) Ae kB T Pt (x ) ∼ tα (16) with α > 0. The logic behind this ansatz is that, instead of Eq. (15) which is obtained from Eq. (14) where the left-hand side is exactly zero, we now look for a time-dependent solution to Eq. (13), where the contribution from the left-hand side is nonzero, yet small. This long-time behavior, when inserted into the Fokker-Planck equation, is a solution to leading-order, with correction terms of order ∂t t −α ∝ t −(1+α ) which are smaller by a factor of 1/t than the leading t −α term. Physically, we √ can expect this solution to be valid only for x ≪ l(t), where l (t ) = 4Dt is the diffusion length-scale of the problem. In the range x ≫ 1, on the other hand, we know that the force is negligible. Hence, in that case, Pt (x ) ≃ √ 1 π Dt exp[− x2 ]. 4Dt (17) 16,17) in the overlap region 1 ≪ x ≪ l(t), we find α = Matching Eqs. (√ 1/2 and A = 1/ π D. A uniform approximation then reads Pt (x ) ≃ √ 1 π Dt e − Vk (xT) − x2 4Dt B e , (18) where we set V (∞ ) = 0. This scaling solution is valid at long times for all x. For a process with a potential of the form in Eq. (12), where the√particle is allowed to cross also to x < 0, the factor √ π Dt → 2 π Dt , and similarly A → A/2 in Eq. (16). In Section 5 we derive this solution for any potential which decays faster than 1/x2 at large distances, using an eigenfunction expansion method that employs the continuous spectrum of the Fokker-Planck operator. This √ method also yields correction terms that vanish faster than 1/ t . We leave the equivalent derivation for the case of potentials that fall off like 1/|x|ζ , where 1 < ζ ≤ 2, out of this manuscript, since here all the results associated with the leading order behavior in time are similar to the ζ > 2 case, but the correction terms are different. Importantly, to obtain Eq. (18), we have assumed that the particles are initially localized, say on x0 or within an interval (0, x0 ), or more correctly the initial density has at least an exponential cutoff. The scale x0 does not alter the long-time solution (to leadingorder approximation, see Section 5), and since the solution forgets its initial conditions, we introduce a thermodynamic notation Pt (x ) ∝ e−βV (x )−ξ x 2 (19) ∞ Fig. 3. (color online) The Gibbs entropy, S (t ) = − 0 Pt (x ) log [Pt (x )]dx, obtained from simulation results of the overdamped Langevin Eq. (1), and the symmetric potential Eq. (12), for 105 particles, at kB T = 0.2 (green squares). The purple circles corresponds to Eq. (21) (but with π → 4π , since the potential is symmetric), where all the observables are measure from the simulation. Here, one should note √ that additional correction terms of order 1/ t might exist, due to contribution from correction terms to the leading-order behavior of Pt (x), discussed in Section 5. Such correction terms will decay in time as fast as V. The blue line corresponds to the asymptotic theory, Eq. (22) (with π → 4π ), which becomes exact when t → ∞. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) with β = 1/kB T and ξ = 1/4Dt. We now consider the entropy of the system. Using the uniform approximation, Eq. (18), we find Sunif (t ) = kB This yields Sunif (t ) =  0 V  T ∞ dxPt (x ) ln √ π Dt + βV (x ) + ξ x2 . + kB ξ x2  − F˜ /T , (20) (21) where the averages are with respect to Eq. (18). The average energy V and β are thermodynamic pairs, and at least suggestively the same is true for x2  and ξ . Later, we will make this analogy to equilibrium statistical mechanics more precise, using an entropy maximization formalism. Here, the free energy is F˜ = −(kB T /2 ) ln(π /4ξ ). It stems from the first term in Eq. (20), which in the usual setting gives the connection between the Helmholtz free energy and the partition function. We now see that, for a fixed observation time and D; Eq. (21) means that 1/T = (∂ Sunif /∂ V  ), in agreement with standard thermodynamics. To leading order, using the fact that the mean-squared displacement is behaving diffusively at long times, namely kB x2 /4Dt = kB /2, we have Sunif (t ) ∼ kB kB − F˜ /T = ln(π eDt ). 2 2 (22) In this limit, the entropy is insensitive to the potential, since the average potential energy approaches zero at increasing times. This occurs simply because the particles increasingly explore the large x region where the potential is flat (zero). Below, we study the average potential energy in detail, but for now it is only important to realize that entropy times T is far larger. Eq. (22) shows that the entropy is increasing with time, which is to be expected, since the packet of particles is spreading out to the medium. Fig. 3 shows the entropy versus time, obtained from simulation results using the two-sided potential, Eq. (12) (green squares), and the corresponding measurement based on Eq. (21) (purple circles). But in the latter, since the potential is symmetric, k Sunif (t ) ∼ 2B ln (4eπ Dt ) (the reason is that at the tails, the density √ is now proportional to exp[−x2 /(4Dt )]/ 4π Dt ). It also shows the asymptotic logarithmic growth at long times (blue line), based on Eq. (22), but with π → 4π . The figure shows that Eq. (21) is a good approximation, but note that additional correction terms of order √ 1/ t might exist, due to contribution from correction terms to the 5 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 leading-order behavior of Pt (x), discussed in Section 5, which will decay in time as fast as V. Even in this case, the asymptotic logarithmic growth at long times, seen in Eq. (22), will remain unchanged. We are now ready to define the non-normalized Boltzmann density more precisely. We define the time dependent function Zt = exp(−F˜ /kB T ), (23) √ √ hence in our case Zt = π Dt (or 4π Dt for two-sided potentials). Inserting Zt into Eq. (18), we find lim Zt Pt (x ) = exp[−V (x )/kB T ]. (24) t→∞ Here we used limt→∞ exp[−x2 /4Dt] = 1 for any finite x, though in a √ finite-time experiment, this identity will be valid for x ≪ l (t ) = 4Dt . At least in principle, with a measurement of the entropy in the long-time limit, which is possible with an ensemble of particles, and using Eqs. (22,23), we can determine Zt . And with that information, we can verify Eq. (24). Note that, clearly, if we know D up front, we can find Zt without measuring entropy at all, but one point of our work is to claim that the thermodynamic formalism may have a more general validity, which is a question worthy of further study. Since Zt is increasing with time, then when it is multiplied by the normalized density Pt (x), this leads to a non-normalizable Boltzmann factor, when t becomes large. Again, mathematically, the non-normalized Boltzmann factor is the infinite-invariant density of the system [2,4]. Of course Eq. (24) holds also for the case when Zt eventually approaches a constant, namely for finite size systems. Fig. 4, shows simulation results for Zt Pt (x ), corresponding to the potentials in Eq. (11) (top figure), and Eq. (12) (bottom). At increasing times, the scaled distributions approach the non-normalizable Boltzmann √ state, Eq. (24), while the finite-time Gaussian tails where x ∼ t , are pushed increasingly outward. The insets show the potentials. Note that we obtain the correction terms to the leading behavior of Pt (x) in Section 5, but for now we leave the rest of the discussion about the corrections to the mean energy and entropy out of the scope of this manuscript. The benefit of that discussion is in the rigorous derivation of the shape of Pt (x) and the effect of the initial condition, which is shown to be negligible (namely it does not alter Eq. (24)). Using Eq. (24), one can determine the shape of the potential field in the system from the position density of the particles. We do not, however, address the question of whether this yields the mechanical or electrostatic force, or an effective force, as clearly the potential of the mean force might itself be temperaturedependent. 4. Entropy extremum principle The structure of the theory suggests that a more general principle is at work. The entropy extremum principle is a natural choice, with the imposition of three added constraints. These are: the normalization condition, a finite mean energy condition (this allows us to treat fluctuations of energy in the canonical ensemble, unlike the microcanonical ensemble, where the energy surface is fixed), and the special feature of our system, which is that the meansquared displacement is diffusive, x2  ∼ 2Dt. This latter constraint is the new ingredient of the theory and this introduces the time dependency. We define a functional of the density Pt (x) at some fixed large t S[Pt (x )] = − kB  ∞ 0 − β kB  Pt (x ) ln Pt (x )dx ∞ 0 V (x )Pt (x )dx − V   Fig. 4. (Color online) Upper panel: The non-normalizable Boltzmann state Eq. (24), corresponding to the Lennard-Jones potential, Eq. (11); Theory in black√line, simulation results for Zt Pt (x ), at various increasing times t, where Zt = π Dt (kB T = 0.0436, γ = 1), appear in colored symbols. The inset shows the potential, where U0 = 10 0 0 and (a, b) = (1, 2 ). At finite times one can see deviations from the asymptotic Boltzmann state, which vanish as t increases. As seen by the yellow line (at e.g., t = 2500), at large x, Pt (x) is Gaussian (since the attractive force is zero). Lower panel: The non-normalizable Boltzmann state corresponding to the two-sided potential Eq. (12), with (kB T, γ ) = (0.8, 1 ) (green curve). Here Zt → 2Zt . Simulation results in colored symbols (and the inset shows the potential). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) − λkB − ξ kb  ∞ Pt (x )dx − 1 0  0 ∞   x2 Pt (x )dx − 2Dt . (25) Here β , λ and ξ are Lagrange multipliers. The first term is the entropy, which in the absence of the constraints implies that all micro-states are equally probable. If we set ξ = 0, and find the extremum, we get the usual Boltzmann-Gibbs theory, however that can be valid only if the potential is binding, which is not the case under study here. Taking the functional derivative, we get   Pt (x ) = N exp −β V (x ) − ξ x2 , (26) where N is the normalization constant. The constraints are N N N    ∞ 0 ∞ 0 ∞ 0 exp −β V (x ) − ξ x2 dx = 1, x2 exp −β V (x ) − ξ x2 dx = 2Dt, and V (x ) exp −β V (x ) − ξ x2 dx = V (x ). (27) These conditions are satisfied if ξ is small and V(x) is asymptotically flat. To see this, we use ξ x2 = y2 , so we can rewrite the first two integrals as: N  ξ  ∞ 0  exp −β V (y/   ξ ) − y2 dy = 1 (28) 6 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 and N  ξ 3  ∞ 0  y2 exp −β V (y/    ξ ) − y2 dy = 2Dt, (29) and since limξ →0 V (y/ ξ ) = 0 we may ignore the potential field in this limit. Solving the Gaussian integrals, we find ξ= 1 , 4Dt N= √ 1 π Dt (30) . Notice that in Eqs. (28,29), we are averaging over observables which are non-integrable with respect to the Boltzmann infinite density, i.e. x0 , x2 (or y0 , y2 ). For the last constraint we use the small parameter ξ to approximate V (x ) exp[−β V (x ) − ξ x2 ] → V (x ) exp[−β V (x )], and since in any reasonable range where the potential is finite, exp[−ξ x2 ] is equal to 1, we get V  = ∞ N 0 V (x ) exp[−β V (x )]dx. To summarize, we find that the extremum principle gives exp −β V (x ) − x2 /4Dt Pt (x ) ≃ . √ π Dt (31) This is the same as the uniform approximation Eq. (18), which was proven valid for a specific stochastic model, i.e. the overdamped Langevin equation. However, the extremum principle suggests that the existence of the non-normalized Boltzmann density is of more general validity, even outside of this particular context. Finally, one could claim that since thermodynamics is a theory which does not depend explicitly on time t, we cannot identify β with the inverse of the temperature. However, at least within our Langevin model, the Einstein relation and our analytical results give both the physical and the mathematical motivation to make this relation. 5. Eigenfunction expansion, and corrections to the uniform approximation We have obtained the leading-order behavior of Pt (x) at long times for asymptotically flat potentials from two different directions, first by using physically inspired guesses for the small and large x regimes and then matching, and secondly via entropy maximization. Here, we re-derive our result a third time, but importantly, now we use an eigenfunction expansion, so as to allow access to the leading-order corrections. In particular, we focus on potentials that fall off faster than 1/x2 at large x. The spectrum of the Fokker-Planck operator is continuous since the system is diffusive in the bulk. For convenience, we consider the case that there is a reflecting wall at x = 0, giving rise to a no-flux boundary condition, Pt′ (0 ) = 0. The final answer works as well for the case that the potential diverges to +∞ as x → 0 from above, so that again the particles are restricted to x > 0. For a δ -function initial distribution, centered around some positive x0 , we show that the initial condition does not affect the asymptotic shape to leading order in time, and we obtain the correction to the leading-order term where it first makes its appearance. Note that the eigenfunction expansion in the case of a two-sided potential follows along the same lines, but given the details of the setup one may need to examine both symmetric and non-symmetric solutions for the eigenmodes. Starting from the Fokker-Planck equation, ∂t Pt (x ) = LˆFP Pt (x ), where erator (32) LˆFP = D[∂x2 + ∂xU ′ (x )] is and U (x ) = Vk(xT) , using B the the Fokker-Planck opunitary transforma- tion [5] Hˆ = eU (x )/2 LˆF P e−U (x )/2 , (x, t ) = eU (x)/2 Pt (x ), we write the following Schrödinger-like equation Hˆ (x ) = −λ (x ). (33) For technical reasons, we imagine an infinite effective-potential wall at large x = L, so that the spectrum is discrete. As the eigenvalues λ are positive definite (λ = 0 is not in the spectrum as it is not normalizable in the L → ∞ limit, and hence the system does not reach equilibrium [5]), we write λ ≡ Dk2 , for some discrete set of k’s. We will eventually take the L → ∞ limit, before we take the large t limit. Then, Pt (x ) = e−U (x )/2+U (x0 )/2  Nk2 k 2 ( x0 ) k (34) (x )e−Dk t . {k} It will prove convenient to set the normalization of k (x) via the condition k (0 ) = e−U (0 )/2 , so we have to incorporate an explicit normalization factor Nk2 . Thus, Eq. (33) translates to ∂2 ∂ x2 k (x ) + U ′′ (x ) 2 k (x ) − U ′2 ( x ) 4 k (x ) = −k2 k ( x ). (35) with boundary conditions k ′ ′ (0 ) = − U (0 ) e−U (0 )/2 . (0 ) = e−U (0)/2 , (36) 2 k The long-time limit is clearly dominated by the small-k modes, so we need to consider only them. We need to treat two regimes separately, first the range x ≪ 1/k, where the right-hand side of Eq. (35) is always small, (denoted region I), and second, for x ≫ 1 (region III). We match the two asymptotic limits in the overlap region 1 ≪ x ≪ 1/k (region II). In Region I, the term −k2 k (x ) is negligible due to the smallness of k. To leading order, then, we have the homogeneous equation, ∂2 ∂ x2 k (x ) + U ′′ (x ) 2 k (x ) − U ′2 ( x ) 4 k ( x ) = 0, (37) with the solution (satisfying the no-flux boundary condi−U (x )/2 , corresponding to tions Eq. (36) at x = 0) h (x ) = e the non-normalizable zero-mode. To next order, we write 2 −U (x )/2 f (x ). Plugging this ansatz into Eq. (35), k (x ) ∼ h (x ) − k e we get −U ′ (x ) f ′ (x ) + f ′′ (x ) = 1. (38) The boundary conditions, Eq. (36), which apply for any k, translate to f (0 ) = f ′ (0 ) = 0 and so a simple calculation yields f (x ) =  x ′ eU (x ) 0  x′ 0 ′′ e−U (x ) dx′′ dx′ . (39) The behavior of f(x) for large x can be analyzed as follows. Define x ′ f1 (x ) ≡ 0 e−U (x ) dx. Then, for large x,  x  ′  e−U (x ) − 1 + 1 dx 0   ∞  ′ =x+ e−U (x ) − 1 dx − f 1 (x ) = 0 ≈ x + ℓ0 ∞ x   ′ e−U (x ) − 1 dx ∞ (40) ′ (e−U (x ) − 1 )dx (note where we have defined the length ℓ0 ≡ 0 that ℓ0 is related to the second virial coefficient, see Section 8). Now, for large x, f (x ) =  x ′ eU (x ) f1 (x ) − x − ℓ0 + x + ℓ0 dx 0  ∞ x2 ′ + ℓ0 x + eU (x ) f1 (x ) − x − ℓ0 dx = 2 0  ∞ ′ eU (x ) f1 (x ) − x − ℓ0 dx − x x2 ≈ + ℓ0 x + A 2 (41) where the constant A with units of length2 is defined as A ≡  ∞ U ( x′ ) e f1 (x ) − x − ℓ0 dx. This behavior of f can be seen to be 0 7 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 5. (color online) Left panel: The potential U (x ) = (1 − 4x2 )/(1 + x6 ), (blue line), and the effective potential in the Schrödinger Eq. (35): = −U ′′ (x )/2 + U ′2 (x )/4 (magenta dashed line). Right panel: The eigenfunction k (x), Eq. (45) (red diamonds), compared with the exact numerical solution of the Schrödinger Eq. (35) (black line), with the potential U (x ) = (1 − 4x2 )/(1 + x6 ) and k = 0.088. At small x, the eigenfunction reflects the shape of the potential, since it is proportional to exp[−U (x )/2] (region I). At large x (region III), the solution is proportional to Ak cos(kx ) + Bk sin(kx ) (see Eq. (43)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) consistent with the differential equation, Eq. (38). In the matching region II, where 1 ≪ x ≪ 1/k, therefore: match ( x ) ≈ 1 − k2 ( x2 /2 + ℓ0 x + A ) (42) (x ) ∼ −k2 (43) In region III, since x ≫ 1, the U′′ (x) and U′ 2 (x) terms are negligible, and therefore, Eq. (35) now reads ∂2 ∂ x2 k k ( x ),  √  √  π Dt Pt (x ) compared with the prediction from Eq. (47), 12 e−U (x ) ( f (x ) + f (x0 ) − 2A + ℓ20 ), (labelled t = ∞ ) for the case U (x ) = (1 − 4x2 )/(1 + x6 ) with x0 = 3, for t = 40, 80, 160, 320. D = 1. For this Fig. 6. The predicted scaled correction Dt e−U (x ) − potential, ℓ0 = 2.1719, A = −5.3464. whose solution is k (x ) ≈ Ak cos(kx ) + Bk sin(kx ). Comparing this solution with Eq. (42) in the matching region, we find that Ak = 1 − k2 A and Bk = −kℓ0 . We see that, if U (x ) = 0, we have l0 = 0, and thus Bk = 0, namely in the force-free case the sin (kx) term is absent, as it should. We are now in a position to calculate the normalization Nk2 = 2 2 ≈ (1 + 2Ak2 − ℓ20 k2 ). L L(A2k + B2k ) (44) It is interesting to note that the presence of the Bk term, induced by the presence of U(x), results in a O(1) leftward shift of ℓ0 in the original pure cos wave of the free particle case, in additional to the small change in normalization. This has a simple physical interpretation, which we will return to below. A uniform approximation, for any x, is seen to be:  −U (x )/2 cos(kx )(1 − k2 ( f˜(x ) + A )) k (x ) ≈ e  −kℓ0 sin(kx ) 1 −U (x ) e 2 6 (45) where we have defined f˜(x ) ≡ f (x ) − x2 /2 − ℓ0 x − A. (46) In the left panel of Fig. 5, we show the potential, U (x ) = (1 − 4x2 )/(1 + x6 ), and the effective potential in the “quantum” problem, namely: −U ′′ (x )/2 + U ′2 (x )/4. In the right panel, we show that k (x) from Eq. (45) matches the exact numerical solution of the Schrödinger Eq. (35), with the above mentioned potential and k = 0.088. We are now in a position to take the L → ∞ limit, wherein the   sum over n transforms into an integral over k, k → dk πL . For finite x, and x0 , then, using Eqs. (34) and (44,45) the longtime density reads Pt (x ) ≈ e−U (x )  ∞ 0 2dk π [1 − k2 ( f (x ) + f (x0 ) − 2A + ℓ20 )] 2 ×e−Dk t   1 e−U (x ) 1− ( f (x ) + f (x0 ) − 2A + ℓ20 ) = √ 2Dt π Dt  π Dt Pt (x ) compared with the ( f (x ) + f (x0 ) − 2A + ℓ0 ), (labelled t = ∞ ) for the Fig. 7. The predicted scaled correction Dt e−U (x ) − prediction from Eq. (47), case U (x ) = (1 − 4x2 )/(1 + x ) with x0 = 1 and x0 = 3, for t = 320. D = 1, ℓ0 = 2.1719, A = −5.3464. giving us the zeroth-order time-dependent Boltzmann-Gibbs factor, with a 1/t correction that grows quadratically in x. We test this prediction  √ in Fig.  6, where we plot the scaled correction Dt eU (x ) − π Dt Pt (x ) and the prediction from Eq. (47), 12 e−U (x ) ( f (x ) + f (x0 ) − 2A + ℓ20 ), for the case U (x ) = (1 − 4x2 )/(1 + x6 ) with x0 = 3. We see that the numerics is converging to the prediction with increasing t, with the size of the correction growing as x2 . To test the dependence on x0 , we plot in Fig. 7, the predicted correction, and the simulation results at t = 320 for the same potential, for both x0 = 1 and x0 = 3. The formula is seen to correctly capture the x0 dependence, with the correction larger in magnitude for larger x0 , √ For large x, of order t , but x0 still O(1), the long-time density is Pt (x ) ≈ (47)  ∞ 0 2dk π   1 − k2 A cos(kx ) − kℓ0 sin(kx ) × 1 − k2 f (x0 ) 1 + 2k2 A − k2 ℓ20 e−Dk 2 t 8 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 More generally, for arbitrary x, we have the uniform solution 2 Px (t ) ≈ e−(x+ℓ0 ) /4Dt e−U (x ) √ π Dt × 1− √ Fig. 8. The √ predicted scaled correction in the outer regime, x ∼ Dt , √ −x2 /4Dt Dt (e − π Dt Pt (x )) compared with the prediction from Eq. (49), 2 ℓ0 x (−A + ℓ20 + f (x0 ))e−x /4Dt , (labelled t = ∞ ) for the case U (x ) = (1 − 4x2 )/(1 + x6 ) 2 with x0 = 1, for t = 40, 160, 640, and 2560. D = 1, ℓ0 = 2.1719. 2Dt − x2 (−A + ℓ20 /2 + f (x0 ) + f˜(x )) , 4 D2 t 2 (50) where f˜(x ) is defined in Eq. (46). In Section 8.1, we use the next-to leading order behavior of Pt (x) in order to obtain a correction term to the leading-order, linear behavior of the mean-squared displacement x2 . Comment on normalization. We can verify that the solution in Eq. (50) is normalized to unity, in the following way: Integrating exp[−U (x ) − x2 /(4Dt )], we find that √  ∞ −U (x )− x2  x2 4Dt dx ≈ π Dt + 0∞ (e−U (x) − 1 )e− 4Dt dx, from which, 0 e 2 at large x, in the longfor potentials that fall-off faster than 1/x√ ∞ 2 time limit we get 0 e−U (x )−x /4Dt dx ≈ π Dt √ + ℓ0 . So, the first ∞ term of 0 Pt (x )dx is approximately 1 + ℓ0 / π Dt . Similarly, the ∞ √ 2 ℓ0 x/2Dt term in Eq. (49), yields − 0 ℓ0 xe−x /(4Dt ) dx/[2 π (Dt )3/2 ]. This cancels out the correction to the normalization from the leading term. 6. Time and ensemble-averages Let us now focus on the limit of long times, where the correction terms to the leading-order behavior of Pt (x) are negligible with respect to the uniform approximation, Eq. (18). To define the long-time limit of averages, we distinguish between two types of observables: integrable and non-integrable observables, with respect to the non-normalized Boltzmann state. We consider first the indicator function O[x(t )] = I (x1 < x(t ) < x2 ), √ Fig. 9. The predicted second order correction in the outer regime, x ∼ Dt ,   scaled √ √ ℓ0 x −x2 /4Dt Dt (e 1 − 2Dt − π Dt Pt (x )) compared with the prediction from Eq. (48), 2 2 2Dt−x −x /4Dt e , (labelled t = ∞ ) for the case U (x ) = (1 − 4x2 )/(1 + x6 ) with x0 = 1, 4D2 t 2 for t = 40, 160, 640, and 2560. D = 1, ℓ0 = 2.1719, A = −5.3464. . e−x = √ 2 2Dt − x2 ℓ x 1− (−A + ℓ20 + f (x0 )) − 0 2Dt 4 D2 t 2 /4Dt π Dt e−x Pt (x ) ≈ √ 2 /4Dt π Dt  1−  ℓ0 x . 2Dt where I (· · · ) = 1 if the condition in the parenthesis is satisfied, and zero otherwise. Along the trajectory x(t) of the particle, this observable, I( · ), switches between values +1 and 0, corresponding to whether the particle is present in the domain (x1 , x2 ) or not. Here, x1 and x2 are the experimentalist’s matter of choice. The ensemble-average of this observable, which in principle can be obtained from a packet of non-interacting particles, at some time t is I (x1 < x(t ) < x2 ) = ∼ (48) Thus, it turns out that in this regime, the leading order correction to Pt comes from a rightward shift in the Gaussian by an amount ℓ0 , due to the shift in the phase of k discussed above. Thus, at large distances, the position of the diffusive √ source is effectively at x = −ℓ0 . As this shift leads to an O(1/ t ) relative change in the solution, if we consider just this leading change, the O(1/t) terms are negligible, and the solution simplifies to (49) Note that ℓ0 may take either positive or negative values hence the sign of the correction term depends on the force field (see Section 8, where we relate ℓ0 to the virial theorem). This prediction is tested in Fig. 8, where we see very good agreement. In Fig. 9, we check the validity of the 1/t relative change via the difference of the simulation to the first order outer solution, Eq. (49), and see that here too the agreement is excellent. (51)  ∞ 0  x2 I (x1 < x < x2 )Pt (x )dx e−βV (x ) /Zt . (52) x1 This result is valid in the limit of long times, when √ x1 and x2 are much smaller than the diffusion length-scale l (t ) = 4Dt , namely we used the approximation exp[−(x2 )2 /4Dt] ≃ 1. We see that, while the Boltzmann factor exp(−β V (x )) is not normalized, it is used to obtain the ensemble averages. In this case the observable is zero at large distances, hence this observable cures the nonintegrability of the infinite density. More generally, for observables integrable with respect to the non-normalized Boltzmann factor we have, using Eq. (18) lim Zt O (x ) = t→∞  ∞ 0 exp[−β V (x )]O (x )dx. (53) Eqs. (52,53) are valid also for the case when the system reaches a steady state, and then Zt is the normalizing partition function; in that case Eq. (52) is simply the probability of finding the particle in thermal equilibrium in the interval (x1 , x2 ). The time that a particle spends in the domain (x1 , x2 ) is called the residence time or the occupation time, and it is denoted tx1 <x<x2 . This variable fluctuates from one trajectory to another, E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 however when the system reaches a steady state (i.e. if the potential is binding), the occupation fraction in the long measurementtime limit clearly satisfies limt→∞ tx1 <x<x2 /t = Prob(x1 < x < x2 ) and the latter is obtained from the Boltzmann measure  x2 x1 lim tx1 <x<x2 /t = t→∞ e −β V ( x ) d x Z (54) . This result can be obtained also from Birkhoff’s ergodic hypothesis in standard Boltzmann-Gibbs theory. What is the corresponding behavior for weakly binding potentials where infinite ergodic theory is relevant? The observable tx1 <x<x2 /t is the finite-time average tx1 <x<x2 = t t 0 I (x1 < x(t ) < x2 )dt . t (55) Let us first consider the ensemble-average of this observable, which is obtained by averaging over an ensemble of paths, each trajectory yielding its own residence time. Here we have t x1 <x<x2 t  =   0t I (x1 < x(t ) < x2 )dt  t . (56) To calculate this value we can switch the order of the ensembleaveraging procedure with the time integration, and use I (x1 < x x(t ) < x2 ) = x 2 Pt (x )dx. Now, we need to perform the time in1 tegration, however considering the long-time limit (and neglecting short-time effects), this calculation is straight forward: using √ Zt = π Dt , Eq. (23), we get t x1 <x<x2 t  1 ∼ t =2  t dt 0  x2 x1  x2 x1 e −β V ( x ) d x Zt exp[−β V (x )]dx Zt . (57) The factor 2 is a consequence of the time integration, since  t −1/2 dt = 2t 1/2 , and note that we may take here the lower limit 0t of the integration to zero, without any affect on the long-time limit. As for the indicator function, now consider the averaged potential energy, with the uniform approximation, Eq. (18): V (x ) ∼ ∞ 0 exp(−β V (x ) − x2 /4Dt )V (x )dx . Zt (58) In the long-time limit, we have lim Zt V (x ) = t→∞  ∞ 0 exp(−β V (x ))V (x )dx, (59) since the potential is zero beyond some length-scale l1 , and exp(−x2 /4Dt ) ∼ 1 for 0 < x < l1 . The total potential energy is decreasing with time (in absolute value, and in contrast with the entropy which is increasing), as particles are escaping the well, traveling to the bulk and exploring the spatial domain where the force is negligible. Here, it is important to note that the process is recurrent, so any particle which escapes the surface to any distance as long as we wish, will eventually return to the regime of nonzero potential with probability one (to the local minimum of the Lennard-Jones potential, for example). This means that if we perform an experiment with N ≫ 1 particles, it is more likely to find them in the medium, beyond l1 , after some finite time. Still, since one always observes the return of the particles, there is always a non-negligible number of them which are residing in the vicinity of the surface (no escape is forever). Now, consider the ensemble-average of the deterministic part of the force field f (x ) = −∂xV (x )  f ( x ) = kB T ∞ 0 ∂x e−V (x)/kB T dx Zt kB T = {exp[−βV (∞ )] − exp[−βV (0 )]}. Zt (60) 9 Since this observable is also integrable with respect to the infinite density, Eq. (24), we get  f (x ) = kb T /Zt or lim √ t→∞ π Dt  f (x ) = kB T , (61) if we consider the case of the one-sided system, and = 0.5kB T in the two-sided case (since Zt → 2Zt ). Notice that this limit does not depend on the specific shape of the potential. For all the integrable observables above (and in fact for any integrable observable), we find a connection between the time and ensemble-average, which is a generalization of the Birkhoff law from standard thermodynamics, namely the doubling effect seen in Eq. (57) is a general feature for this class of observables. Consider an observable O[x(t )] which is integrable with respect to the non-normalized Boltzmann factor, then the ensemble-average of the time-average is O[x(t )] = 2O (x ) (62) where O (x ) = ∞ 0 O (x ) exp[−β V (x )]dx . Zt (63) The factor 2 is a consequence of the diffusive nature of the process, which leads to the integration over the time-dependent partition function, hence this doubling effect might be widely observed. The numerical results which support Eq. (62), were presented in our previous work, Ref. [2], where we showed that the simulations agreed with the theory. Below, in Section 9.2, we show numerical evidence for a variation of Eq. (62), which is valid also in dimensions d ≥ 1. 7. Fluctuations of the time-averages The time-average O[x(t )] in the time-independent canonical setting is equal to the ensemble-average, in the long-time limit (ergodicity). In our case, the time-averages fluctuate between different trajectories, which is a common theme in single-molecule experiments, and here we explore the fluctuations. To start, we again consider the indicator function, O[x(t )] = I (x1 < x(t ) < x2 ), defined in Eq. (51). However, our results are far more general than that. As we will show, the fluctuations of time-averages of observables integrable with respect to the non-normalized Boltzmann density follow a universal law, in the spirit of the Aaronson-Darling-Kac theorem [4]. For simplicity, let us consider x1 = 0. The process I(x1 < x(t) < x2 ), starting inside the region (0, x2 ), is switching randomly between two states, with sojourn times in the interval close to the surface denoted τ in , when I (x1 < x(t ) < x2 ) = 1, and τ out , when I (x1 < x(t ) < x2 ) = 0. The first time interval in the domain (0, x2 ) is denoted τ1in , and the rest follow, so the sojourn times in the two states are given by the sequence   τ1in , τ1out , τ2in , τ2out , · · · . (64) These times can be treated as mutually independent, identically distributed random variables, since temporal correlations in the Langevin Eq. (1) decrease exponentially fast in time. We denote the probability density functions of out and in sojourn times with ψ out/in (τ ), respectively. Importantly, in the long-time limit; ψout (τ ) ∝ τ −3/2 . This well-known result is related of course to the flatness of the potential field at large x. In this regime, the process x(t) is controlled by diffusion and while it is recurrent, so the density ψ out (τ ) is normalized, the average sojourn time in the out state diverges. This absence of a typical timescale, together with the diverging partition function, are precisely the reasons for the failure of the standard (Birkhoff) ergodic theory, and the emergence of Boltzmann-like infinite ergodic theory. The sojourn times in the in state are thinly distributed and, importantly, the moments 10 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 of τ in are finite, which is clearly the case since the interval (0, x2 ) is of finite length. Let n be the number of switching events from the in to the out states. For a fixed measurement time t, this number is random. We claim that the distribution of n is determined by the statistics of the out times, when t is large. Roughly speaking, the in sojourn times are very short when compared to the out times, since those have an infinite mean. This means that in the time interval (0, t), we will typically observe an out sojourn time of the order of magnitude of the measurement time t, and the size of this largest interval controls the number of out-to-in transitions (if the largest out time is very long, then n is small, compared with a realization with a shorter maximal out time). Note that the time-average of I(x1 < x(t) < x2 ) is equal to the sum of the in times, divided by the measurement time; I = n in i=1 τi /t, where we assume, without any loss of generality, that at time t the process is in state out. The mean I was already obtained rigorously from the non-normalized density in Eq. (57), but with the notations of renewal theory we have I ≃ τ in n/t and since n ∼ t1/2 we have I ∝ t −1/2 , as we found earlier. The behavior n∝t1/2 is well known in renewal theory, and with a handwaving argument we note that the effective average sojourn time t is τeff  = 0 τ ψout (τ )dτ ∝ t 1/2 , and hence n ∼ t/τ eff ∝t1/2 . Let us now consider the second moment of the time-average. Since, for any i = j, (τiin )2  = (τ jin )2  and τiin τ jin  = τiin 2 , we argue that 2 I  = =  in 2 i=1 i t2  n τ   n(τ in )2  + n(n − 1 )τ in 2 t2 . (65) Now we are ready to get to the main point of this section. Considering the variance of the time-average  ( I ) 2  − I 2 = n 2  − n 2 t2 (tr )2  − tr 2 τ in 2 + t2 n  t2 = (τ in )2  − τ in 2 ,    (66) (67) Var(τ in ) we see that the second term is negligible, compared with the first, since from renewal theory we know that n ∼ t1/2 . This means that the fluctuations of τ in are irrelevant, and there is a single important timescale describing the process, which is the average sojourn time τ in . We now normalize the variance using I ∼ nτ in /t and find 2 n 2  − n 2 I  − I 2 . → n 2 I 2 (68) This analysis can be continued to higher order moments and it yields I I  → n n  ≡ η. (69) This means that the residence time in (x1 , x2 ), divided by the mean residence time, is equal in a statistical sense to number of the into-out transitions over their mean. The probability density function of 0 < η is known from renewal theory [41], and since ψout (τ ) ∝ τ −3/2 we find PDF(η ) = 2 π 2 e − η /π . (70) This has, by definition, unit mean. Naively, the reader might be tempted to believe that this result is related to the Gaussian central limit theorem. However, this is not the case, since for sojourntime distributions with other fat tails we will get a form very different from Gaussian [2,16,42]. In-fact, PDF(η) is a special case of a more general density function known as the Mittag-Leffler distribution, which in turn is related to Lévy statistics. Note, that the Aaronson-Darlin-Kac theorem [4] predicts that the distribution of the time average of a process with an infinite measure, will be given by the Mittag-Leffler distribution in the form of Mα (η ) =  Ŵ 1/α (1+α ) Ŵ 1/α (1+α ) , with lα ,0 (# ) being the one-sided Lévy l αη1+1/α α ,0 η1/α density (defined as the inverse-Laplace transform of exp(−uα ), from u → #), see e.g., [22,43,44]. The exponent α , we argue, is determined by the first return probability ψout (τ ) ∝ τ −1−α , and in our case, α = 1/2 means that M1/2 (η ) is equal to Eq. (70). Other values of α are found in the case of diffusion in logarithmic potentials, as explained in Ref. [2] (see also [5]), which are out of the scope of this manuscript, but in that case a derivation of the Mittag-Leffler distribution can also be made following the same lines as presented below. A hand-waving argument for Eq. (70), works as follows: Consider n independent, identically distributed random variables τ1out , . . . τnout , which correspond in our physical model to the times in the state out. According to the Lévy central limit theorem [45], the probability density function (PDF) of these times is the one-sided Lévy density with index 1/2 l1/2,0 ( τ ) = 2 1 √ π τ −3/2 exp − 1 . 4τ  (71) Here, like ψ out (τ ), l1/2,0 (τ ) ∝ τ −3/2 and this fat-tailed behavior allows us to consider a specific choice of the out times distribution, in the sense that asymptotically the results are not sensitive to the short τ behavior of ψ out (τ ). We use dimensionless units, and since eventually we consider the dimensionless variable η, this is not ∞ a problem. The Laplace transform LT [ f (τ )] = 0 f (τ ) exp(−τ u )dτ 1 / 2 of Eq. (71) is exp(−u ). Now, consider the random variable y = n 2 i=1 τi /n . The PDF of y is also the one-sided stable law Eq. (71), since it is easy to check that exp(−uy ) = exp(−u1/2 ). We are interested in the probability distribution of n, and we fix the mea surement time t to be the sum of the sojourn times t = ni=1 τi . 2 Hence y = t/n , and PDF(n ) = PDF(y )| l1/2,0 t n2 dy |= dn  2  2t 1 n | 3 | = √ exp − . 4t n πt (72) Thus the density of n is half a Gaussian. From here, we find n = √ 2 t/π , and switching to the random variable η = n/n we get Eq. (70). Throughout this derivation, we treat n as a continuous variable, which makes sense in the long-time limit, and can be justified using well-known rigorous results. We note that, mathematically, the number of switching events n is formally infinite. This is related to the fact that the Langevin trajectories are continuous, hence once we have one transition, we experience many of them. This is not a major problem since we actually considered the scaled random variable η = n/n which has a unit mean. To put it differently, since n ∼ t1/2 in the longtime limit, we consider a scaled variable which is perfectly well behaved. From the measurement point of view, we sample the trajectory with a finite rate, so n is always finite, and this is also true in simulations, where we use discrete steps in space and time (in the limit of large n, the results will not be sensitive to the sampling rate and the discretization). Inspired by infinite ergodic theory, we claim that the MittagLeffler distribution of time-averages is a far more general result. For example, consider the time-average of the potential energy V . Also here the observable V[x(t)] is switching between long periods where it is nearly zero (when the particle is far in the bulk), to relatively short bursts when this observable is nonzero, when 11 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 11. The value of Zt xF (x ), obtained from simulation results based on numerical integration over the Langevin Eq. (1), where F (x ) = −V ′ (x ), and V (x ) = (−25x2 + x4 /2 ) exp(−x2 /20 )/125 (green circles), approaches at increasing √ times to the theoretical value given by Eq. (75). Here Zt = 2 π Dt , since the potential is symmetric. γ = 1, and D = kB T = 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the potential energy, using various numerical simulations, whose results are presented in our previous paper, Ref. [2]. 8. Virial theorem Fig. 10. Upper panel: Simulated time series of the potential energy of a particle in the Lennard-Jones potential Eq. (11). It exhibits long periods during which |V[x(t)]| < ǫ , where ǫ > 0 can be as small as we wish, and short periods where the absolute value of the energy is high. Small absolute values of the potential energy corresponds to events where the particle strays far away from the potential minimum (the particle being in the range x ≫ 1). Lower panel: The probability distribution g(τ ) of the durations (τ s) of the time periods where the absolute value of the energy is lower than ǫ = 10−5 (here, the total measurement time was t = 2 × 106 ). log [g(τ )] versus log (τ ), obtained from the simulation (magenta circles), is matched with a linear fit (orange line) with a slope of ≈ −3/2, which implies that g(τ ) has a power-law shape similar to free Brownian motion. the particle is close to the surface. Again the statistics of the number of times that the particle visits the domain where V(x) is nonnegligible, is similar to that of n, and its statistics is controlled by the first-passage probability density function from the bulk to the vicinity of the wall. Again, the latter is the fat-tailed density with the familiar τ −3/2 law that we have just mentioned above. So we have O[x(t )] O[x(t )] → η, (73) and Eq. (70) still holds. Note that this yields a complete description of the problem in the sense that O is calculated in principle with the non-normalized Boltzmann density, and we assume that the observable is integrable with respect to this state. In the upper panel of Fig. 10, we see a simulated sample of the time series of the potential energy, of a particle in a Lennard-Jones potential (the details of the simulation are similar to Fig. 4). The time series exhibits long periods where |V(x)| < ǫ , and ǫ is some lower cutoff that can be as small as we wish, and short periods of high energy (in absolute value). In the lower panel of Fig. 10, we see that the probability distribution g(τ ) of the durations (τ s) of the events where the energy is low, which correspond to events where the particle has strayed far away from the potential minimum, has a power-law shape g(τ ) ∝ τ −3/2 , at large τ (as seen from the fitting function, in an orange line), like in free Brownian motion, as expected. Hence the process is recurrent, but the meanreturn time is infinite. In this example, we used ǫ = 10−5 . We verified the Mittag-Leffler distribution of several observables, including The virial theorem addresses the mean of the observable xF(x), where F (x ) = −V ′ (x ). binding potentials, treated with standard thermodynamics, yield xF (x ) = −kB T . In our case [2], using the non-normalized Boltzmann state we find, by integration by parts xF (x ) = kB T Zt =−  kB T Zt ∞ x 0  ∞ 0  ∂  −V (x)/kB T e − 1 dx ∂x   e−V (x )/kB T − 1 dx, (74) where we used our convention exp[−V (∞ )/kB T ] = 1. Now, using ∞ ℓ0 = 0 [exp(−V (x )/kB T ) − 1]dx (introduced in Section 5), we get lim Zt xF (x ) → −kB T ℓ0 . t→∞ (75) The ratio ℓ0 /Zt distinguishes Eq. (75) from the standard thermal virial theorem, where the ratio −xF (x )/kB T at equilibrium is unity. Note that ℓ0 = −2B2 , (76) where B2 is called the second virial coefficient [1]. For twosided, symmetric potentials xF (x ) = −2kB T ℓ0 /Zt (but as mentioned, also Zt → 2Zt ). The constant l0 points to a surprising link between the virial theorem and the corrections to the uniform approximation, studied in Section 5 (particularly, Eq. (49) for large x), which means that by measuring the shape of the tails of the diffusing particle packet, one can, at least in principle, obtain knowledge about the force in the system, even though in the large x region it is effectively zero. Interestingly, notice that ℓ0 can change sign, for various potentials, which is also very different from standard thermodynamics. Fig. 11 shows the approach of the simulated value of Zt xF (x ) where F (x ) = −V ′ (x ), and V (x ) = (−25x2 + x4 /2 ) exp(−x2 /20 )/125 at increasing times (green circles), to the theoretical limit, Eq. (75) (red line) with kB T = γ = 0.5, where the value is negative. This result was obtained from the overdamped Langevin Eq. (1). In Fig. 12, we show the various values of l0 , obtained for the potential U (x, U0 ) = U0 (x/5 )4 − (x/5 )2 exp (−x/5 )2 , for various values of the amplitude U0 (blue line), at fixed temperature: kB T = 1, as well as for a fixed U0 = 2, and varying T (green line). In both those cases we used γ = 1. 12 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 12. (Color online). Various values of l0 , obtained for the potential U (x, U0 ) = U0 (−x/5 )4 − (x/5 )2 exp (−x/5 )2 , which change sign for various values of the amplitude U0 (blue line), at fixed temperature: kB T = 1, as well as for a fixed U0 = 2, and varying T (green line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) kB T = γ = 1, which confirm the validity of the correction term to the leading-order, linear, behavior of the second moment in time. The figure also shows an additional constant coming from higherorder correction terms, which was obtained numerically. Note that in the case of a two-sided system, there might be additional√ correction terms to the mean-squared displacement, of order t , if the initial position of the particle is not located at the origin. The reason is that, here, the correction terms to Pt (x) might differ from Section 5. To understand the connection between Eq. (77) and the virial theorem, we need to extend our analysis to the phase space and consider both the particle’s position x, and it’s velocity v. Namely, in what follows, instead of Eq. (1) we now use the underdamped Langevin equation, with zero-mean, white Gaussian noise; m∂t v = F (x ) − γ v + η (t ). In this process, which also obeys the fluctuationdissipation relation, γ > 0, and we include the acceleration term according to Newton’s second law, where m is the particle’s mass. The analysis below will yield the same results regardless if V(x) is a one-sided or two-sided potential, given that the process starts at x(t = 0 ) = 0. Consider the identity m∂t xv = mv2  + mx∂t v. (78) m∂t xv = kB T + xF (x ) − γ xv, (79) Since the velocity is thermalised mv2  = kB T , and since xη (t ) = 0 we get where clearly xv = ∂t x2  2 . (80) Furthermore, we assume that in the long-time limit Fig. 13. (Color online) Numerical results (blue line), corresponding to overdamped Langevin dynamics with the same potential used for all the figures in Section 5, with kB T = γ = 1), which confirm the validity of the correction term to the leadingorder, linear, behavior of the second moment in time. The fitting curve is shown in √ orange-dashed line, and it corresponds to the t term in Eq. (77), plus an additional constant coming from higher-order correction terms, which was obtained numerically. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 8.1. A Correction Term for the mean-squared displacement, and the underdamped langevin process The result of the previous subsection provides us with a nice example that demonstrates that the next-to-leading order correction term, derived for the uniform approximation of Pt (x) in Section 5, has also some relevance to thermodynamics. This link is made by looking at the correction term to the linear behavior, in time, of the mean-squared displacement x2  of the system, which we now derive using an elementary calculation, and then explain it in terms of the virial theorem. Using Eq. (49), we get x2  ≈  ∞ x2 π Dt 0 − √  ∞ 0 2 e−x /(4Dt ) dx 2   x2 e−x /(4Dt ) ℓ0 x dx √ 2Dt π Dt = 2Dt − 4ℓ0  Dt/π , x2  ∼ 2Dt + Dβ t β + . . . , (81) γ xv = kB T + xF (x ). (82) where the second term is small compared with the first. It is then clear that |m∂t xv| ≪ |γ xv| and hence For equilibrium situations, i.e. for binding potentials like the Harmonic oscillator, xv = 0 since the marginal position density is described by the Boltzmann distribution, and Maxwell-distribution for the velocities. This means that in thermal equilibrium, the velocities are not correlated with the spatial position of the particle, since in the single particle Hamiltonian the kinetic energy is separated from the potential energy. For the case under study here, the correlation xv is not strictly zero. Coming back to Eq. (82), we see that the term xF(x) decays like t −1/2 , since Zt ∝ t 1/2 . Using the leading order term xv = ∂t x2 /2 ∼ D, from Eq. (82) we get in the very long-time limit kB T − γ D = 0, recovering the Einstein relation. Hence we need to consider the sub-leading terms. It is easy to see that since xF (x ) ∝ t −1/2 we must have β = 1/2. Using Eqs. (80,81,82) we find xF (x ) = D1/2 γ , 4t 1/2 (83) and from Eq. (75) (77) where we used the fact that at large x, V(x) ≈ 0. In the above derivation, we considered the integration limits to be zero and infinity, although Eq. (49) is exact only in the large x limit, since the contribution to the mean-squared displacement from the small x regime is negligibly small at the long-time limit. Fig. 13 shows numerical results corresponding to overdamped Langevin dynamics with the same potential used for all the figures in Section 5, with D1/2 √ 4 Dℓ0 . =− √ π (84) It follows that D1/2 is determined by the potential energy via the length-scale ℓ0 , and for a given D it is independent of the mass of the particle. Using Eq. (84) in Eq. (81), we therefore recover the first and second leading order terms in x2 , obtained in Eq. (77). Fig. 14 shows the approach of the mean xF(x) obtained from simulation results, using the underdamped Langevin equation with the potential Eq. (12), to the asymptotic value Eq. (75). 13 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 yields ∂ ∂ 1 − d ∂2 P (r ) = D + 2 Pt (r ), ∂t t ∂r r ∂r (87) ∞ where 0 Pt (r )dr = 1. The solution to this equation, for various boundary conditions, is found e.g. in [47–49]. In two dimensions, for example, starting from a ring-shaped initial distribution P0 (r, φ ) = δ (r − r0 )(0 < φ < 2π )/(2π r0 ), with a reflecting boundary condition at r = 0 (see details in [47]), at time t we find  Pt (r ) ≈ (r /2Dt ) exp −(r 2 + r02 )/4Dt I0 (r r0 /2Dt ), where Iν ( · ) is the Bessel function with index ν (and here ν = 0). In d-dimensions, starting from a uniform probability distribution on a d-dimensional sphere of radius r0 , and a reflecting boundary condition at r = 0 [47]: Fig. 14. (Color online) The approach of Zt xF (x ) obtained from simulation results corresponding to the underdamped Langevin process with the two-sided potential Eq. (12), with kB T = 0.8 and γ = 0.5, at increasing times (blue circles) to the asymptotic theoretical value Eq. (75) (green line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) We remark that the observable xv is non-integrable with respect to the non-normalized Boltzmann-Gibbs state. In phase space, we speculate that the non-normalized state is exp[−H/kB T ]/Zt where the Hamiltonian H = p2 /2m + V (x ) and p = mv as usual. The only change  here is that Zt is now √ π Dt 2π mkB T , where the factor 2π mkB T stems from the Maxwell distribution. This expression describes the bulk fluctuations of the packet of particles in phase space, while in the far tails the correlation between x and v builds up. The full analysis of the phase-space infinite-density remains out of the scope of this manuscript. 9. Non-normalizable boltzmann-Gibbs states in d-dimensions So far, we have treated only one-dimensional processes. However, as we mentioned in the introduction, the issue of the nonnormalizability of the Boltzmann factor raised by Fermi was in the context of three-dimensional motion, under the influence of a Coulomb potential [40]. We now show that the non-normalizable Boltzmann-Gibbs state is found also in d-dimensions, when the external potential is isotropic, and it decays at least as rapidly as 1/r, at large distances. One should keep in mind here that in 2dimensions, in the absence of any potential field, Polya’s theorem states that a Brownian particle still returns to its origin with probability 1 (and the mean first return time is infinite), but in any dimension d > 2, this is no longer the case, if the system size is infinite. This holds also in the presence of Coulomb-like potentials. Still, as we now show, the Boltzmann infinite density is valid. We begin our analysis in the absence of any force. The radial motion of a Brownian particle in d-dimensions, in the space defined by the orthogonal directions χ1 , χ2 , . . . χd , is described by Bessel process [46] d−1 √ r˙ = D + 2DŴ(t ), r where r =  (85) χ12 + χ22 + . . . χd2 . Accordingly, the radial Fokker-Planck equation describing the expansion of the probability density Wt (r) in d dimensions is [46] ∂t Wt (r ) = D ((d − 1 )/r )∂r + ∂r2 Wt (r ), ∞ where 0 Wt (r )1r d−1 dr = 1 and c(d) is a constant which rises from the integration over all the angular degrees of freedom of the ddimensional Laplacian. Here we assumed that the initial distribution of the particles was also isotropic around the origin. Substituting Wt (r ) = Pt (r )r 1−d /c (d ), (86) Pt (r ) = r0 2Dt r r0 d/2  exp − r 2 + r02 4Dt  Id/2−1  r r0 . 2Dt (88) Below, we add to the system an asymptotically flat potential field, with a repelling part on the origin. In that case, the right-hand side of Eq. (88) will describe the shape of the density in the large r regime, where V(r) ∼ 0. One can easily verify, for example, that when d = 1 and r0 ≪ t1/2 , but r ∼ O (t 1/2 ), Pt (r) in Eq. (88) is ∼ √ 1 exp(−r 2 /(4Dt )) (as is well known [47]). π Dt 9.1. The infinite-invariant density Consider an isotropic, radially dependent potential V(r), such that V (0 ) = ∞ and V (∞ ) = 0, which falls-off at least as rapidly as 1/r at large distances (in the same spirit as the potentials we investigated in the unidimensional case). Now, the radial dynamics is described by r˙ = D (d − 1 ) r − V ′ (r ) γ + √ 2DŴ(t ). (89) As in the unidimensional Langevin equation, Eq. (1), here D = kB T /γ and kT , γ > 0. The corresponding radial Fokker-Planck equation is ∂ V ′ (r ) ∂2 d−1 ∂ ∂ Wt (r ) = D + r 1−d r d−1 + 2 Wt (r ). ∂t r ∂r ∂r kB T ∂r (90) Here, and in what follows, we assume that the initial particle distribution is narrow, and isotropic. After repeating the substitution from Eq. (86), this yields ∂ ∂ 1 − d ∂ V ′ (r ) ∂ 2 Pt (r ) = D + + 2 Pt (r ). ∂t ∂r r ∂ r kB T ∂r (91) To solve Eq. (91) to leading order, in the long-time limit, we use the ansatz Pt (r ) ≈ r0 2Dt r r0 d/2 exp − r 2 +r02 4Dt  Id/2−1  rr0  2Dt I (r ), where I (r ) is some general function of r. Note that this approach is similar to that employed in the unidimensional case, in Section 3. Plugging this ansatz in Eq. (91), we obtain the uniform approximation Pt (r ) ≈ r0 2Dt r r0 d/2 e− r 2 +r 2 0 4Dt Id/2−1  r r0 − V (r ) e kB T , 2Dt (92) for long t. From the uniform approximation, Eq. (92), using the asymptotic shape  rr  of the Bessel  function in the  limit t → ∞, since Id/2−1 2Dt0 ≈ (rr0 )d/2−1 / (4Dt )d/2−1 Ŵ(d/2 ) , and exp[−(r 2 + r02 )/(4Dt )] → 1, we find lim Zt r 1−d Pt (r ) → exp (−V (r )/kB T ), t→∞ (93) where Zt = 2d/2−1 Ŵ(d/2 )(2Dt )d/2 . Importantly, from this relation we again see that, in the long time limit, the non-normalizable solution is independent of r0 . In one dimension, from this result we recover Eq. (18). 14 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 15 shows excellent agreement between the simulation results of a two-dimensional Langevin process, and the theory corresponding to Eqs. (92,93). The simulation method used an EulerMayurama integration √ scheme over the two Langevin √equations x˙ = −V ′ (r ) cos(θ ) + 2DŴ1 (t ) and y˙ = −V ′ (r ) sin(θ ) + 2DŴ2 (t ), where V (r ) = (r /5 )4 − (r /5 )2 exp[−(r/5 )2 ]. (94) V (r ) = (2/r )12 − 0.5/r. (95) Ŵ 1 (t) and Ŵ 2 (t) represent two independent, zero-mean and δ correlated Gaussian white noise terms. At t = 0, the particles were uniformly distributed around a ring of radius r0 = 5. Fig. 16 shows simulation results of a three-dimensional Langevin process, with the Coulomb-type potential Here r is defined in spherical coordinates. The repulsive part of the potentials, which falls-off as rapidly as 1/r12 , was added to regularize the interactions near the origin in the simulation (this is technically simpler to realize numerically than putting a hard spherical wall around some r ≪ 1). This model mimics a field created by a finite-sized charge, which repels the observed particle at short distances. The figure shows excellent match between the simulated radial PDF Pt (r), multiplied by Zt (defined in Eq. (93)) at times t = 50 0 0, 50 0 0 0, 10 0 0 0 0 (green stars, blue circles and red diamonds, respectively), and the corresponding uniform approximation, Eq. (92) with d = 3 (magenta, purple and brown lines, respectively). At increasing times, the simulation results approach the non-normalizable Boltzmann state exp(−V (r )/kB T ) (black line), via Eq. (93), as expected, confirming the existence of the infiniteinvariant density also in three dimensions. Here, D = kB T = 0.05. In this system the minimum of the field is on rmin ≈ 2.84, hence in Fig. 16 we see a peak at this value. Also notice that the depth of the well is V (rmin ) ≈ −0.16, hence V (rmin )/kB T ≈ −3.22. Thus we are dealing here with a trap that is not too deep, this allows the escape of the particles on a finite observation time. For a deeper well, we will need to wait for even longer observation time to find the infinite density. Fig. 15. (Color online) Simulation results for a two-dimensional Langevin process, with V(r) in Eq. (94) (γ = 1, and  kB T = D = 0.5), compared with theory. Colored symbols: (Zt /r )Pt (r ), where r = x2 + y2 , Zt = 2Dt (see Eq. (93)) and Pt (r) is the radial distribution, obtained from the simulation at times t = 50 0 0, 250 0 0, 46322 (blue circles, orange squares and magenta circles, respectively). The corresponding theoretical curves obtained from Eq. (92) with d = 2 (multiplied by 2Dt/r) are the red dash-dot line and a green dashed line and cyan dashed-line. The nonnormalizable Boltzmann factor exp(−V (r )/kB T ), appears in solid black line. Here, the initial distribution of the 106 particles was uniform around a ring of radius r0 = 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) now relate this weighted time average to an ensemble average, performed at time t. Clearly, Zt O (r )t = Zt Therefore, lim t→∞ 9.2. Ergodicity of time-weighted observables, in d-dimensions As mentioned above, Brownian motion in d-dimensions is nonrecurrent. In this section, we propose a new approach for evaluating time and ensemble averages of observables integrable with respect to the infinite-density, Eq. (93), and show that this method leads to a Birkhoff-like equality between the two means, which is valid in any dimension d ≥ 1. Let O (r ) be an integrable observable, with respect to r d−1 exp (−V (r )/kB T ), in Eq. (93), e.g., O (r ) = (ra < r (t ) < rb ). Here, (ra < r (t ) < rb ) = 1 while the particle’s trajectory passes inside the d-dimensional shell with inner radius ra > 0, and outer radius rb < ∞, and zero otherwise. We define the ensemble mean of the weighted time-average as   Zt O [r (t )] ≡ = 1 t  t 0   1 t 0 t Zt˜ O [r (t˜)]dt˜ Zt˜ O (r )t˜ dt˜ = 1 t  0 t dt˜Zt˜  0 ∞ O (r )Pt˜ (r )dr. (96)   ∞ 1 t 1 dt˜Zt˜ O (r )r d−1 e−V (r )/kB T dr, t 0 Zt˜ 0  ∞ = O (r )r d−1 e−V (r )/kB T dr. Zt O [r (t )] → ∞ 0 O (r )Pt (r )dr → O (r )nB .  Zt O [r (t )] Zt O (r )t → 1. (98) (99) Eq. (99) is a generalized form of Birkhoff’s theorem. Note that in one dimension, this equation constitutes an alternative to Eq. (62), which was derived for non-weighted time-averages. Here, like in Section 6, both the ensemble-average, and the weighted time-average are estimated in the long-time limit from the nonnormalizable (d-dimensional) Boltzmann-Gibbs state, Eq. (93), hence Eq. (99) constitutes further extension of infinite-ergodic theory (in the sense of the ratio between time and ensemble means). Fig. 17 shows the ratio between the weighted time average Eq. (96), measured in a two-dimensional Langevin simulation, and the ensemble average with respect to the non-normalizable Boltzmann state nB of the indicator function (2.5 < r(t) < 5). This ratio converges to unity at increasing times, validating Eq. (99). The details of the simulation are similar to Fig. 15. Fig. 18 shows a similar ratio as in Fig. 17, but this time it is measured in a threedimensional Langevin simulation, for (2.8 < r(t) < 3) (blue line). The details of this simulation are similar to Fig. 16, and the ratio again converges to unity at increasing times. 10. Summary of our main results Using Eq. (93), in the limit t → ∞, Eq. (96) leads to     (97) 0  We denote O (r )nB ≡ 0∞ O (r )rd−1 e−V (r )/kB Tdr, and then Eq. (97) means that when t → ∞, Zt O [r (t )] → O (r )nB . We In this manuscript, which extends our work in Ref. [2], we have shown that the spatial shape of a diffusing particle packet, inside an asymptotically flat potential, converges at long times to a nonnormalizable Boltzmann state in any dimension, d ≥ 1. This state, which is treated mathematically as an infinite-invariant density [4], takes the place of the standard Boltzmann distribution, which gives the equilibrium state in systems with strong confinement, in E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Fig. 16. (Color online) Simulation results of a three-dimensional Langevin process, with the Coulomb-type potential Eq. (95). The figure shows excellent match between the simulated radial PDF Pt (r) (here we used spherical coordinates), multiplied by Zt /r 2 (and Zt is defined below Eq. (93)) at times t = 50 0 0, 50 0 0 0, 10 0 0 0 0 (green stars, blue circles and red diamonds, respectively), and the corresponding uniform approximation, Eq. (92) with d = 3 (magenta, purple and brown lines, respectively). At increasing times, the simulation results approach to the nonnormalizable Boltzmann state exp(−V (r )/kB T ) (black line), via Eq. (93), as expected, confirming the existence of the infinite-invariant density also in three dimensions. Here, D = kB T = 0.05, and we used 107 particles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 15 the sense that we can use it to obtain time and ensemble averages of integrable observables. We have mainly focused on onedimensional systems, which obey the Aaronson-Darlin-Kac theorem, and here we also showed that the non-normalizable Boltzmann state gives the entropy-energy relation, and the virial theorem. We studied the emergence of the latter in detail in one dimension, and it would be interesting to study it further also in higher dimensions in a future work. We have obtained the non-normalizable Boltzmann state in one dimension using three different techniques: via physical scaling assumptions (Section 3), using the entropy maximization principle (Section 4), and via a rigorous eigenfunction expansion method (Section 5). The last of these also yielded terms which describe the sub-leading order behavior of the probability density function, which decay more rapidly with time. Though the analysis based on an eigenfunction expansion in d > 1 dimension is left for future work, here we showed that by attaching a proper weight function to integrable observables, the ratio between weighted time averages and ensemble averages converges to unity (see Eq. (99)). The distribution of the weighted time average is an open question for future research. The main results of this manuscript are the uniformly valid approximation for the one-dimensional probability density Pt (x) including the first-order correction, Eq. (50); the relationship between the mean-squared position and the virial theorem, expressed in Eq. (77) and Eqs. (81 and 84); the leading order probability density in arbitrary dimensions, Eq. (93); and lastly, the ergodicity of time-weighted observables, expressed in Eq. (99). 11. Discussion Fig. 17. (Color online) The ratio between the weighted time average Eq. (96), measured in a two-dimensional Langevin simulation, and the ensemble average with respect to the non-normalizable Boltzmann state nB of the indicator function (2.5 < r(t) < 5), converges to unity at increasing times (simulation results in magenta line, the green line on 1 serves as a guide to the eye). This validates Eq. (99). The details of the simulation are similar to those in Fig. 15 (so nB ≈ 11.7076). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 18. (Color online) The ratio between the weighted time average Eq. (96), measured in a three-dimensional Langevin simulation, and the ensemble average with respect to the non-normalizable Boltzmann state nB of the indicator function (2.8 < r(t) < 3), converges to unity at increasing times (simulation results in blue line, the green line on 1 serves as a guide to the eye). This validates Eq. (99). The details of the simulation are similar to Fig. 16 (so nB ≈ 41.7144), except that here we used 105 particles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Infinite ergodic theory can be applied to many thermodynamic systems, as the main condition is that the fluctuation-dissipation theorem holds. One must distinguish however between recurrent, and non-recurrent processes, since only in the latter the Aaronson-Darlin-Kac theorem holds. Physically, the key point is that we can identify easily important observables that are integrable with respect to the infinite density, and the fluctuationdissipation theorem guarantees that this infinite density is the non-normalizable Boltzmann factor. Ryabov, et al. [9,10] considered a different, though related, setup in one dimension, with an unstable potential that does not allow for the return of the particle to its starting point. Also here the partition function diverges, but again a certain aspect of Boltzmann equilibrium remains. Hence, our work suggests that further aspects of ergodicity should be studied in this setup, perhaps in the spirit of the evaluation of weightedtime-averaged observables, and it also indicates that future studies in that direction could be interesting in other classes of potentials as well. Our work also encourages the investigation of these problems in non-Markovian settings, and in the presence of manybody interactions, and since we have assumed that the fluctuationdissipation relation holds, it will be interesting to examine if this assumption can be relaxed and infinite-ergodic theory can be studied also e.g. in the framework of active particles, as in [50,51]. While in this manuscript we considered the fluctuations of time averages, and in particular the fluctuations of energy, the whole framework of stochastic thermodynamics could in principle be investigated. For example, it would be interesting to explore the fluctuations of the rate of entropy production, and the work and heat exchange between the particles and the heat bath. It should be noted however that our theory gives rise to both extensions of stochastic thermodynamics, for systems with a non-normalizable Boltzmann-Gibbs state, but also the connection between fluctuations (diffusivity) and thermodynamics. This is seen in the virial correction to the diffusion law (Section 8). In the current theory, one cannot separate diffusion from non-normalized Boltzmann- 16 E. Aghion, D.A. Kessler and E. Barkai / Chaos, Solitons and Fractals 138 (2020) 109890 Gibbs states, as was demonstrated in the extremum principle studied in Section 4. Credit Author Statement All the authors contributed to theory in the manuscript, and various calculations. DAK and EA performed also simulations. All the authors also contributed equally to the writing of the manuscript together. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Acknowledgement The support of Israel Science Foundation grant 1898/17 is acknowledged. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.chaos.2020.109890 References [1] Chandler D. Introduction to modern statistical mechanics. Introduction to Modern Statistical Mechanics, by David Chandler, pp 288 Foreword by David Chandler Oxford University Press, Sep 1987 ISBN-10: 0195042778 ISBN-13: 9780195042771 1987:288. [2] Aghion E, Kessler DA, Barkai E. From non-normalizable Boltzmann-gibbs statistics to infinite-ergodic theory. Phys. Rev. Lett. 2019;122(1):010601. [3] Darling D, Kac M. On occupation times for markoff processes. Trans. Am. Math. Soc. 1957;84(2):444–58. [4] Aaronson J. “An introduction to infinite ergodic Theory”. American Mathematical Soc.; 1997. [5] Dechant A, Lutz E, Barkai E, Kessler DA. Solution of the fokker-planck equation with a logarithmic potential. J. Stat. Phys. 2011;145(6):1524–45. [6] Bouin E, Dolbeault J, Schmeiser C. Diffusion and kinetic transport with very weak confinement. Kinetic & Related Models 2020;13(2):345. [7] Wang X, Deng W, Chen Y. Ergodic properties of heterogeneous diffusion processes in a potential well. J. Chem. Phys. 2019;150(16):164121. [8] Sivan M, Farago O. Probability distribution of brownian motion in periodic potentials. Physical Review E 2018;98(5):052117. [9] Šiler M, Ornigotti L, Brzobohatỳ O, Jákl P, Ryabov A, Holubec V, et al. Diffusing up the hill: dynamics and equipartition in highly unstable systems. Phys. Rev. Lett. 2018;121(23):230601. [10] Ryabov A, Holubec V, Berestneva E. Living on the edge of instability. J. Stat. Mech.: Theory Exp. 2019:084014. [11] Zweimüller R. Surrey notes on infinite ergodic theory. Lecture notes, Surrey Univ 2009. [12] Thaler M. Infinite ergodic theory. course notes “from the dynamic odyssey”, cirm 2001. 2001. [13] Akimoto T, Miyaguchi T. Role of infinite invariant measure in deterministic subdiffusion. Physical Review E 2010;82(3):030102. [14] Korabel N, Barkai E. Infinite invariant density determines statistics of time averages for weak chaos. Phys. Rev. Lett. 2012;108(6):060604. [15] Klages R. Weak chaos, infinite ergodic theory, and anomalous dynamics. In: From Hamiltonian Chaos to Complex Systems. Springer; 2013. p. 3–42. [16] Akimoto T, Shinkai S, Aizawa Y. Distributional behavior of time averages of non- l∧ 1 l1 observables in one-dimensional intermittent maps with infinite invariant measures. J. Stat. Phys. 2015;158(2):476–93. [17] Meyer P, Kantz H. Infinite invariant densities due to intermittency in a nonlinear oscillator. Physical Review E 2017;96(2):022217. [18] Leibovich N, Barkai E. Infinite ergodic theory for heterogeneous diffusion processes. Physical Review E 2019;99(4):042138. [19] Akimoto T, Barkai E, Radons G. Infinite invariant density in a semi-markov process with continuous state variables. arXiv preprint arXiv:190810501 2019. [20] Zhou T, Xu P, Deng W. Continuous time random walks and l\’{e} vy walks with stochastic resetting. arXiv preprint arXiv:190907213 2019. [21] Sato Y, Klages R. Anomalous diffusion in random dynamical systems. Phys. Rev. Lett. 2019;122(17):174101. [22] Radice M, Onofri M, Artuso R, Pozzoli G. Statistics of occupation times and connection to local properties of nonhomogeneous random walks. Physical Review E 2020;101(4):042103. [23] Kessler DA, Barkai E. Infinite covariant density for diffusion in logarithmic potentials and optical lattices. Phys. Rev. Lett. 2010;105(12):120602. [24] Lutz E, Renzoni F. Beyond Boltzmann-gibbs statistical mechanics in optical lattices. Nat. Phys. 2013;9(10):615–19. [25] Rebenshtok A, Denisov S, Hänggi P, Barkai E. Non-normalizable densities in strong anomalous diffusion: beyond the central limit theorem. Phys. Rev. Lett. 2014;112(11):110601. [26] Rebenshtok A, Denisov S, Hänggi P, Barkai E. Infinite densities for Lévy walks. Phys Rev E 2014;90:062135. [27] Holz PC, Dechant A, Lutz E. Infinite density for cold atoms in shallow optical lattices. EPL (Europhysics Letters) 2015;109(2):23001. [28] Aghion E, Kessler DA, Barkai E. Large fluctuations for spatial diffusion of cold atoms. Phys. Rev. Lett. 2017;118(26):260601. [29] Wang W, Vezzani A, Burioni R, Barkai E. Transport in disordered systems: the single big jump approach. Physical Review Research 2019;1(3):033172. [30] Vezzani A, Barkai E, Burioni R. Single-big-jump principle in physical modeling. Physical Review E 2019;100(1):012108. [31] Pavani SRP, Thompson MA, Biteen JS, Lord SJ, Liu N, Twieg RJ, et al. Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function. Proceedings of the National Academy of Sciences 2009;106(9):2995–9. [32] Chechkin AV, Zaid IM, Lomholt MA, Sokolov IM, Metzler R. Bulk-mediated diffusion on a planar surface: full solution. Physical Review E 2012;86(4):041101. [33] Metzler R, Jeon J-H, Cherstvy AG, Barkai E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. PCCP 2014;16(44):24128–64. [34] Campagnola G, Nepal K, Schroder BW, Peersen OB, Krapf D. Superdiffusive motion of membrane-targeting c2 domains. Sci. Rep. 2015;5:17721. [35] Krapf D, Campagnola G, Nepal K, Peersen OB. Strange kinetics of bulk-mediated diffusion on lipid bilayers. PCCP 2016;18(18):12633–41. [36] Wang D, Wu H, Schwartz DK. Three-dimensional tracking of interfacial hopping diffusion. Physical Review E 2017;119(26):268001. [37] Ashkin A, Dziedzic JM, Bjorkholm J, Chu S. Observation of a single-beam gradient force optical trap for dielectric particles. Opt. Lett. 1986;11(5):288–90. [38] Grier DG, Roichman Y. Holographic optical trapping. Appl. Opt. 2006;45(5):880–7. [39] Drobczynski S, Ślezak J. Time-series methods in analysis of the optical tweezers recordings. Appl Opt 2015;54(23):7106–14. doi:10.1364/AO.54.007106. [40] Fermi E. Über die wahrscheinlichkeit der quantenzustände. Zeitschrift für Physik 1924;26(1):54–6. [41] Godreche C, Luck J. Statistics of the occupation time of renewal processes. J. Stat. Phys. 2001;104(3):489–524. [42] He Y, Burov S, Metzler R, Barkai E. Random time-scale invariant diffusion and transport coefficients. Phys. Rev. Lett. 2008;101(5):058101. [43] Korabel N, Barkai E. Pesin-type identity for intermittent dynamics with a zero lyaponov exponent. Phys. Rev. Lett. 2009;102(5):050601. [44] Akimoto T, Barkai E. Aging generates regular motions in weakly chaotic systems. Physical Review E 2013;87(3):032915. [45] Metzler R, Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 20 0 0;339(1):1–77. [46] Kubo R, Toda M, Hashitsume N. “Statistical physics II: nonequilibrium statistical mechanics”, 31. Springer Science & Business Media; 2012. [47] Bray AJ. Random walks in logarithmic and power-law potentials, nonuniversal persistence, and vortex dynamics in the two-dimensional xy model. Physical Review E 20 0 0;62(1):103. [48] Martin E, Behn U, Germano G. First-passage and first-exit times of a Bessel-like stochastic process. Physical Review E 2011;83(5):051115. [49] Medalion S, Aghion E, Meirovitch H, Barkai E, Kessler DA. Size distribution of ring polymers. Sci. Rep. 2016;6:27661. [50] Romanczuk P, Bär M, Ebeling W, Lindner B, Schimansky-Geier L. Active brownian particles. The European Physical Journal Special Topics 2012;202(1):1–162. [51] Hoell C, Löwen H, Menzel AM. Multi-species dynamical density functional theory for microswimmers: derivation, orientational ordering, trapping potentials, and shear cells. J. Chem. Phys. 2019;151(6):064902.