key: cord-0118630-o57i2dfm
authors: Horii, Hiroshi; Lefevere, Raphael; Itami, Masato; Nemoto, Takahiro
title: Anomalous fluctuations of renewal-reward processes with heavy-tailed distributions
date: 2022-03-16
journal: nan
DOI: nan
sha: 9a4e8b99775dbe70c9413fceb1cb6f79332ad371
doc_id: 118630
cord_uid: o57i2dfm

For renewal-reward processes with a power-law decaying waiting time distribution, anomalously large probabilities are assigned to atypical values of the asymptotic processes. Previous works have reveals that this anomalous scaling causes a singularity in the corresponding large deviation function. In order to further understand this problem, we study in this article the scaling of variance in several renewal-reward processes: counting processes with two different power-law decaying waiting time distributions and a Knudsen gas (a heat conduction model). Through analytical and numerical analyses of these models, we find that the variances show an anomalous scaling when the exponent of the power law is -3. For a counting process with the power-law exponent smaller than -3, this anomalous scaling does not take place: this indicates that the processes only fluctuate around the expectation with an error that is compatible with a standard large deviation scaling. In this case, we argue that anomalous scaling appears in higher order cumulants. Finally, many-body particles interacting through soft-core interactions with the boundary conditions employed in the Knudsen gas are studied using numerical simulations. We observe that the variance scaling becomes normal even though the power-law exponent in the boundary conditions is -3.

A renewal-reward process, a generalisation of continuous time Markov processes, is one of the simplest stochastic processes that can describe random sequences with memory effects [1] [2] [3] . In contrast to its the Markov counterpart, in renewal-reward processes, the waiting time to move from one state to the next one can be distributed by a non-exponential function. The process can thus describe a broad spectrum of phenomena in physics [4] and other fields, including a melt up of the stock market [5, 6] and a super spreader in epidemics [7, 8] , where memory effects are known to be important.

When the waiting time distribution has a power law, the dynamics show a slow convergence to its stationary states due to its heavy tail. For example, the probability that the state of the system always stays in the initial state during the dynamics remains nonnegligible in the large time limit [9] . This anomalous behaviour can be characterised using a large deviation principle (LDP) [10, 11] . LDP states that the logarithmic probability of a time-averaged quantity is proportional with the averaging time (with a negative proportional constant), except for the trivial probability where the time averaged quantity takes its expectation. In renewal reward processes with power-law waiting time distributions, this proportional constant, known as a rate function or large deviation function (LDF), can take the value 0 not only for the expectation but also for a certain range of the values [9, 12, 13] . This indicates that these events are more likely to occur than in standard systems. We call this range of LDF taking the value 0 the affine part.

The affine part tells us that these rare events occur more likely than usual, but does not tell us how likely they do. To solve this problem, finite-time analyses of the LDP are necessary. One such attempt could be a so-called strong LDP, where the next order corrections of the logarithmic probability from the LDP are computed [14] . However, at present, it is not clear how this general theory can be extended to the case with the affine part. In [15] , Tsirelson studied a renewal-reward process with general waiting time distributions and derived the next order correction to the LDP. But he used a condition in which an affine part can not be present. Recently, in [16] , the authors studied finite-time corrections of the moment generating function under the condition that the affine part appears (Theorem 2.1). Yet they did not succeed to translate it to the correction term of the LDP.

In this article, instead of focusing on the probability of rare events, we focus on the variance of the time-averaged quantities. The variance can tell us directly how much the averaged quantities fluctuate. If one considers an exponential function for the waiting time distribution, the variance of the time-averaged quantity decreases proportionally to the inverse of the averaging time because this corresponds to the case of a process having a short memory. This indicates that the averaged value mostly falls in the range around the expectation with an error that is proportional with the inverse square root of the averaging time. In the presence of the affine part when heavy-tailed distributions are used for the waiting times, we identify, in this article, a condition under which this scaling of the variance changes. This is consistent with the fact that heavytailed distributions introduce memory effects. Interestingly, not every power-law decaying distribution will result in this scaling modification of the variance: We show that for distributions whose density decay faster than 1/t 3 , the variance keeps its normal scaling. In that case, we expect that the scaling of higher order cumulants are affected, as discussed at the end of this article.

This article is organised as follows. Two models defined using renewal-reward processes are considered in this article: a counting process and a single particle model of heat conduction. These models are studied in Section II (counting process) and in Section III (heat conduction). Each section is organised with (i) a model introduction, (ii) introduction of renewal equations (a key tool to study the asymptotics of moments), (iii) analyses on the first moment, (iv) analyses on the second moment and variance, and (v) numerical studies. In Section IV, we discuss the scaling in higher order cumulants and how the scaling will change in more general heat-conduction systems. In particular, we observe that when several particles are present and interact through soft-core interactions, the time-average of physical quantities recover a "normal" behaviour. This seems to indicate that interactions break the strong memory effects that are present in the system of non-interacting particles.

A renewal-reward process is a model to describe events that occur sequentially. For a given event, the next event occurs after a random waiting time (also called a renewal time or arrival time). The waiting times are independent-and-identically distributed random positive variables (τ k ) k∈N with a probability density p. For this density, we consider the inverse Rayleigh distribution

and the Pareto distribution

with α = 3, both of which do not have a finite second moment, i.e., E[τ 2 ] = ∞. The main quantity of interest in this section is the number of events that have occured up to time t > 0. This is the counting

where S k = τ 1 + . . . + τ k . We denote its q-th order moment by m q (t):

Note that with respect to [17] , we consider the case where the expectation of the waiting time is finite and the renewal theorem [1] implies that the counting process N t behaves as N t ∼ t/E[τ ] for t → ∞. We study the fluctuations around that behaviour.

To analyse the asymptotics of m q (t), we rely on renewal equations: a powerful tool to analyse renewalreward processes. From a straightforward computation, one can establish the following renewal equation for m 1 (t) [2] :

where 

we then derive, from the equation (5),

where we have usedp(s) = sF (s). Similarly, one can also derive a renewal equation for m 2 (t),

from which the Laplace transform of m 2 (t) is obtained asm 2 (s) =m 1 (s)(1 + 2sm 1 (s)).

Moreover, a renewal equation for the momentgenerating function can be derived.

(See Appendix A). From the equation, we derive the Laplace transform of m q (t) as

C. Convergence of the first moment When a waiting-time density p has a finite mean E[τ ] = µ and a finite variance σ 2 , Feller has proven (Chapter 11, section3, theorem 1) [3] that

This result can be easily derived by using the following expansion:

Indeed, by inserting it into (7), we get

which leads to

A rigolous justification to derive (14) from (13) is based on the Tauberian theorem [3] . See Appendix B for more details. From this argument, we can see that the condition E[τ 2 ] = ∞ is necessary for m 1 (t) to have an anomalous scaling. For this reason, we study in this section the two waiting-time distributions behaving at infinity like 1/t 3 . Let us first consider the case of the inverse Rayleigh distribution. Let

We then have for its cumulative distribution function,F

from (7). We then expand φ(s) in s:

leading tõ

By using the Tauberian theorem (Appendix B), we obtain

for large t. This is to be compared to (11) : we see that the convergence is slower in our case. We can repeat the same analysis in the case of the Pareto distribution. The cumulative distribution is derived as (7) and again look at the expansion around s of m and get

We thus obtain the following behaviour for m(t) for large t:

We then study the large time behaviour of the variance

In the case that a waiting time density p has a finite mean E[τ ] = µ and a finite variance σ 2 , we obtain from (9) and (13) 

which yields

with the aid of the Tauberian theorem (B1). c 2 (t) is finally obtained as

Let us now consider the case of the inverse Rayleigh distribution. Inserting the expression (19) form 1 (s) in (9), we obtain,

and then

Therefore

Proceeding in the same way for the Pareto distribution, we obtain in that case

We perform numerical simulations of the counting process N t to illustrate the accuracy of (20), (23) , (30) and (31) . First, m 1 (t) − t 2/(βπ) (resp. m 1 (t) − t) computed from the numerical simulations is plotted as an orange line in Fig.1 (a) (resp. Fig.1 (b)) for the inverse Rayleigh (resp. Pareto) waiting time distribution. According to (20) and (23) , these lines are equivalent to ln(t)/π + o(ln(t)) and ln(t) + o(ln(t)). Assuming that these o(ln(t)) terms are constant over time when t is large, we next plot ln(t)/π+const. (Fig.1(a) ) and ln(t)+const. (Fig.1(b) ) in the same figures.

We then plot (m 2 (t) − m 1 (t) 2 )/t computed from the same numerical simulations in Fig.1 (c,d) for the inverse Rayleigh ( Fig.1(c) ) and the Pareto ( Fig.1(d) ) waiting time distributions. Reference lines 2 Fig.1(c) ) and 2 ln(t) + const. (Fig.1(d) ) are also plotted in the same figures. In these four figures, we observe good agreements between the slopes of the reference lines and the results of numerical simulations in semi-log scale. This demonstrates the validity of (20), (23), (30) and (31) .

Our aim in this section is to show that the slow convergence of the renewal function of processes having density ∼ 1/t 3 as t → ∞ also holds for physical observables in a Knudsen gas [18] . For this, let us consider the model of a single particle bouncing back between two thermal walls.

We consider a particle in a one-dimensional box that has two different temperatures at both ends. The confined tracer moves freely in the box of size 1 and is reflected at the end of the box with a random speed v distributed according to the following Rayleigh distribution:

and v 0 the initial position and velocity of the particle and σ 0 = v 0 /|v 0 |. We denote the initial condition by θ, i.e., θ = (x 0 , v 0 ). The first time that the particle hits a wall is given by S θ,0 = ( 1 2 (σ 0 + 1) − x 0 )/v 0 , and the subsequent hitting times are given by

where v k is a random variable distributed according to a law q βσ k and σ k = (−1) k σ 0 . This may be rewritten as

with the sequence of independent waiting times (τ k ) k∈N distributed with the inverse Rayleigh distribution p β k (τ ) defined as (1) . The energy exchanged between the two walls during a time interval [0, t] is defined as

if t ≥ S θ 0 and J θ (t) = 0 otherwise, where N t is the counting process (3), We denote by m θ,q (t) the q-th moment of J θ (t):

A generalisation to the system with an arbitrary box size L is straightforward. Indeed, denoting by S L θ,k the corresponding hitting times with the boundaries, it is easy to see that S L θ,k = LS 1 θ,k . This indicates that N L t = N 1 t/L where N L t denotes the counting process corresponding to the hitting times S L θ,k . For the energy current in a box of size L, we also have that β + β − Fig. 2 . Schematic figure to explain the setup of the 1 particle model. When the particle moves to the right (resp. left) wall, σ k = 1 (resp. −1)

In the following we perform all computations with the case L = 1 and then obtain the result for an arbitrary L > 0 by using this scaling relation.

For simplicity, we consider only the following two types of initial conditions:

with v 0 < 0 and

with v 0 > 0, i.e., the cases of a particle just before hitting the left wall (temperature β + ) and of a particle just before hitting the right wall (inverse temperature β − ). As the particle immediately hits each wall when the process starts, the value of the initial velocity v 0 is unimportant. We thus denote by + the initial condition θ + and by − the initial condition θ − .

Dynamics with these two initial conditions are related via the renewal property:

This means that the process conditioned by the firstwaiting time (the left-hand side) is equal to the other process with some increments (the right-hand side). By integrating (40) with respect to the inverse Rayleigh waiting time density (1), we obtain the following coupled renewal-reward equations for the currents

In order to derive the speed of convergence of the current, we perform a Laplace transform of (41) and (42):

whereH ± (s) is the Laplace transform of H ± (t) = 1 2t 2 e − β ± 2t 2 andF β± (s) is the Laplace transform of the cumulative inverse Rayleigh distribution. By substituting (43) into (44), we then obtain an equation for m +,1 (s) as

which leads to

where κ is the conductivity given by

Using again the Tauberian theorem for Laplace transform (Appendix B), we finally get,

Note that the asymptotic form of the average current m +,1 (t) and m −,1 (t) have opposite signs, but this is because the definition of the current includes (−1) ±1 term: these two expressions are physical equivalent. For the average current in a box of size L, we get

(49)

We next discuss the large time asymptotics of the variance of the current. The renewal property for the second moment of the current is expressed by

Then, integrating (50) with respect to the inverse Rayleigh waiting time density (1) and using the relation 

for t > 0, the renewal equations for the second moment of the current are derived as

In order to derive the large time asymptotic of the second moment of the current, we perform Laplace transform of (52) and (53),

(55) We solve these linear equations form −,2 (s) and m +,2 (s)) by usingm ±,1 (s) obtained in section III B and the following expansions ofL ± (s) andg ± (s)

RecallingF β± (s) = β

the Laplace transform of the second moment of the current is derived as

From the Tauberian theorem (Appendix B), we finally arrive at the asymptotic form of the variance

This result agree with our previous work [12] . As for the variance of the current in a box of size L, we get:

(61)

A similar formulation can be applied to study the convergence of the time-averaged kinetic energy defined as

The expected values of the energy are denoted by

As in the section III B, we construct two equations in Laplace spacẽ

Here, h(t) and g(t) are given by

Thus,m + E (s) is calculated as

Proceeding in tha same way as for the current, we can expand the functions involved for small s and obtaiñ

As with the derivation of (48), the large time asymptotics of m ± E (t) are derived as

We numerically simulate the one-particle model to check the validity of (48) and (60). We estimate m +,1 (t) and m +,2 (t) from the numerical simulations, and plot m +,1 (t)−κ(1/β + −1/β − ) and m ±,2 (t)/t 2 − (m ±,1 (t)) 2 /t 2 in Fig. 3 (a,b) . In the same figures, we also plot (κ 2 /2)(β + + β − ) (1/β + − 1/β − ) (ln t)/t + const. and κ 3 (β + + β − ) (1/β + − 1/β − ) 2 + const. as blue dashed lines. We observe that the slopes of the orange lines in semi-log scale asymptotically converge to those of blue dashed lines. This demonstrates (48) and (60).

In the first part of this article, we studied a counting process N t with two heavy-tail waiting time distributions: the Pareto distribution with α = 3 and the inverse Rayleigh distribution. These two waiting time distributions have an asymptotic form 1/τ 3 when the waiting time τ is large, implying that the variance of the waiting time E[τ 2 ] diverges. Because of this divergence, we discussed that the scaled variance c 2 (t)t of the counting process N t also diverges in the large t limit. We indeed derived that it is asymptotically proportional with ln(t), diverging as t → ∞.

A natural question would be, can we get a similar result with a waiting time distribution that has an asymptotic form 1/τ α with α > 3? As demonstrated in Appendix A, one can formulate a general framework, for the Pareto distribution, to derive analytical expressions of the Laplace transform of E[N k t ] (k = 1, 2, 3, ...) for any α. As an example, we computed the first, second and third moments for α = 4, from which we show the third cumulant of N t /t has an asymptotic form ln(t)/t 2 when t is large. This indicates that the third-order cumulant multiplied by 

m+,1(t) − κ(1/β+ − 1/β−) (a) and m±,2(t)/t 2 − (m±,1(t)) 2 /t 2 (b) obtained from numerical simulations (with 10 8 samples) are plotted as orange lines. β+ = 1, β− = 2. Blue dashed lines are (κ 2 /2)(β+ + β−) (1/β+ − 1/β−) (ln t)/t + const. and κ 3 (β+ + β−) (1/β+ − 1/β−) 2 + const. The slopes of the numerical-simulation results in semi-log scale converge to those of the dashed reference lines, showing the validity of (48) and (60). t 2 is asymptotically proportional with ln(t), which is also diverging in the large t limit.

For the counting process with a general fat tail waiting time distribution (that has a power law decay as t → ∞), an existence of the affine part in the scaled cumulant generating function (sCGF) G(s) = lim t→∞ (1/t) ln E[e sNt ] has been proven [16] . When the sCGF is analytic, it can be expanded using scaled cumulantsc i (i = 1, 2, ...) as G(s) = ∞ i=1 (c i /i! )s i by definition, wherec i is defined as lim t→∞ c i t i−1 with the i-th order cumulant c i of N t /t. In the presence of the affine part, sCGF is not analytic, implying that some scaled cumulants lim t→∞ c i t i−1 diverge. Based on the observation above, we conjecture that the k-th order scaled cumulant converges when k < α − 2. When k = α − 2, the k-th order cumulant c i increases proportionally with ln(t)/t i−1 , resulting in ln(t) divergence of lim t→∞ c i t i−1 . It is an interesting future work to study this conjecture.

In the second part of this article, we studied a particle confined in the two walls in different temperatures, and observed that the scaled variance diverges proportionally with ln(t). Here we discuss if we can observe the same divergence in many-body particles confined in the walls.

One-dimensional hard-core interacting particles exchange their velocities when they collide. The dynamics of these particles can thus be exactly mapped to the dynamics of non-interacting many-body particles. Let J N ,L,D ∞ (t) and J N ,L,D 0 (t) be the energy currents of N hard-core interacting and noninteracting particles of diameter D confined in a onedimensional box of size L, respectively. Then, we get

Var

where m L ±,1 (t) and Var(J L ± (t)/t) are given in Eqs. (49) and (61), respectively. This implies that the logarithmic divergence of the scaled variance should be observed in hard-core interacting systems. In soft-core interacting systems, on the other hand, the same mapping cannot be used. This is because of the collisions involving more than two particles, where the exchange rule of velocities no longer holds. To demonstrate this insight, we have performed simulations of hard-core and soft-core interacting particles. The details of the simulations are explained in Appendix C, and the results are shown in Fig. 4 , where J MD (t) is the total energy transferred to the colder wall from time 0 to t, and k is a parameter corresponding to the softness of particles. Note that k = ∞ corresponds to the case of the hard-core interacting system. We observed that the ln(t) divergence disappears as soon as particles start to interact via soft-core interactions. It is an interesting future problem to develop a framework to quantitatively understand the disappearance of the divergence in soft-core particles.

Finally, we list related studies. Studying a variance in a process that is defined with power-law decaying distribution is not something new. In [17] , several anomalous diffusion models were studied using continuous-time random walk, and revealed anomalous scaling in their diffusion coefficients. These anomalous scalings were argued to be universally observed in transports in random media [4] . One of the authors also studied a single big jump principle, which states that the sum of random variables can be approximated by their maximum when the probability distribution of the variables has a power-law [19] .

Singularities of large deviation functions of timecumulative quantities are also known as dynamical phase transitions, and have been studied in many physical models, such as glass formers [20] [21] [22] [23] [24] [25] [26] , lattice gas models [27] [28] [29] [30] [31] [32] [33] , diffusive hydrodynamic equations [34] [35] [36] , and high-dimensional chaotic dynamics [37] [38] [39] and active matters [40] [41] [42] [43] . Finitesize scalings of the large deviation functions have been performed in several works (see an interesting recent work [44] for example), but variance scalings have not been intensively studied in this field yet. 

we calculate the inverse Laplace transform ofm k (s) as m 1 (t) ∼ 2t + 1, (A16) m 2 (t) ∼ 4t 2 + 10t − 16 ln(t), (A17) m 3 (t) ∼ 8t 3 + 48t 2 + 206t − 144t ln(t) (A18) as t → ∞. The second cumulants c 2 defined as (24) and the third cumulant c 3 (defined as the third cumulant of N t /t 3 ) are then given by

(A20) particles is given by

where Θ is the Heaviside step function, and k is a parameter corresponding to the softness of particles. The boundary condition is the same as explained in Sec. III A. Using the second-order symplectic integrator, we numerically solved the equations of mo-tion for the particles, and calculated E[J MD (t)] and Var(J MD (t)) for various values of k. In the simulation, we set the parameter values as N = 3, L = 5, m = D = 1, β + = 1/3, and β − = 1. The timediscretization step-size was set to 0.01. For the case of the hard-core interacting system (denoted by k = ∞), we performed event-driven simulations in which two particles instantaneously exchange velocities when they come into contact. The boundary condition and the parameter values were the same as for the soft-core particle system.

Probability and random processes

An introduction to probability theory and its applications

Packets of diffusing particles exhibit universal exponential tails. Physical review letters

Financial risk and heavy tails

Portfolio value-at-risk with heavytailed risk factors

Evidence that coronavirus superspreading is fat-tailed. Proceedings of the National Academy of

Scale-free networks and sexually transmitted diseases: a description of observed patterns of sexual contacts in Britain and Zimbabwe

Hot scatterers and tracers for the transfer of heat in collisional dynamics

The large deviation approach to statistical mechanics

Large deviations of the current in stochastic collisional dynamics

Macroscopic fluctuation theory of aerogel dynamics

A strong large deviation theorem

From uniform renewal theorem to uniform large and moderate deviations for renewalreward processes

Large time asymptotic of heavy tailed renewal processes

Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking

Model of nonequilibrium ensemble: Knudsen gas

Single-big-jump principle in physical modeling

Dynamic order-disorder in atomistic models of structural glass formers

First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories

Large deviations and ensembles of trajectories in stochastic models

Dynamic transition in an atomic glass former: A molecular-dynamics evidence

Theory of amorphous ices

Finite-size scaling of a first-order dynamical phase transition: Adaptive population dynamics and an effective model

First-order phase transition in a model glass former: Coupling of local structure and dynamics

Distribution of current in nonequilibrium diffusive systems and phase transitions

Universal cumulants of the current in diffusive systems on a ring

Long range correlations and phase transitions in non-equilibrium diffusive systems

Finite size scaling of the dynamical freeenergy in a kinetically constrained model

Dynamical symmetry breaking and phase transitions in driven diffusive systems

Geometrical interpretation of dynamical phase transitions in boundary-driven systems

Universality in dynamical phase transitions of diffusive systems

Current fluctuations in stochastic lattice gases. Physical review letters

Spontaneous symmetry breaking at the fluctuating level

Structure of the optimal path to a fluctuation

Probing rare physical trajectories with Lyapunov weighted dynamics

Large deviations of Lyapunov exponents

Stochastic averaging, large deviations and random transitions for the dynamics of 2D and geostrophic turbulent vortices

Dynamic phase transitions in simple driven kinetic networks

Large fluctuations and dynamic phase transition in a system of selfpropelled particles

Phase separation and large deviations of lattice active matter

Optimizing active work: Dynamical phase transitions, collective motion, and jamming

Finite Time Large Deviations via Matrix Product States

The authors thank S. Sasa for useful comments. Here, we derive the k-th moment of a counting process N t (with a waiting time density p) by using a renewal equation. The moment-generating function M h (t) is defined by