S De Nigris

A Hastir

R Lambiotte

Burstiness and fractional diffusion on complex networks

Many dynamical processes on real world networks display complex temporal patterns as, for instance, a fat-tailed distribution of inter-events times, leading to heterogeneous waiting times between events. In this work, we focus on distributions whose average inter-event time diverges, and study its impact on the dynamics of random walkers on networks. The process can naturally be described, in the long time limit, in terms of Riemann-Liouville fractional derivatives. We show that all the dynamical modes possess, in the asymptotic regime, the same power law relaxation, which implies that the dynamics does not exhibit time-scale separation between modes, and that no mode can be neglected versus another one, even for long times. Our results are then confirmed by numerical simulations.

INTRODUCTION

Linear diffusive processes are often a good entry point for the analysis of the complex dynamics of networked systems. In spite of their apparent simplicity, their flexibility allows for the modelling of a diversity of processes, including information and epidemic spreading [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF]. As an additional assumption, it is also often assumed that the diffusion takes place on a static topology. The dynamical properties of the system are then entirely determined by the spectral properties of a matrix associated to the underlying graph, for instance its Laplacian (for random walks) or its adjacency matrix (for epidemic spreading). In this work, the process can also be used to explore the structure of networks from a dynamical point of view, as exploited in several algorithms for, e.g. community detection or node centrality [START_REF] Lambiotte | Random Walks, Markov Processes and the Multiscale Modular Organization of Complex Networks[END_REF].

Despite its many successes, this simplified model of real-world systems is also known to present several limitations [START_REF] Holme | Temporal networks[END_REF]. In particular, the implicit assumption that events take place at a constant rate, as the realization of a stationary Poisson process, has been contradicted in a variety of systems, in terms of temporal correlations [START_REF] Rosvall | Memory in network flows and its effects on spreading dynamics and community detection[END_REF][START_REF] Scholtes | Causality-driven slow-down and speedup of diffusion in non-Markovian temporal networks[END_REF], non-stationarity [START_REF] Rocha | Bursts of vertex activation and epidemics in evolving networks[END_REF][START_REF] Horváth | Spreading dynamics on networks: the role of burstiness, topology and nonstationarity[END_REF][START_REF] Malmgren | A Poissonian explanation for heavy tails in e-mail communication[END_REF] but also the presence of fattailed distributions of inter-event times [START_REF] Barabási | The origin of bursts and heavy tails in human dynamics[END_REF][START_REF] Eckmann | Entropy of dialogues creates coherent structures in e-mail traffic[END_REF][START_REF] Rocha | Information dynamics shape the sexual networks of Internet-mediated prostitution[END_REF][START_REF] Starnini | Random walks on temporal networks[END_REF][START_REF] Karsai | Small but slow world: How network topology and burstiness slow down spreading[END_REF]. The main purpose of this work is to focus on this latter aspect, and to investigate the impact of the shape of the corresponding distributions on diffusive processes. This question has a long tradition in statistical physics, in situations when the underlying network is a lattice [START_REF] Klafter | First Steps in Random Walks: From Tools to Applications[END_REF]. In recent years, researchers have also been interested in random structures [START_REF] Rocha | Bursts of vertex activation and epidemics in evolving networks[END_REF][START_REF] Vazquez | Impact of non-Poissonian activity patterns on spreading processes[END_REF][START_REF] Iribarren | Impact of human activity patterns on the dynamics of information diffusion[END_REF][START_REF] Jo | Analytically solvable model of spreading dynamics with non-Poissonian[END_REF][START_REF] Min | Spreading dynamics following bursty human activity patterns[END_REF], seen as models of complex networks. It is only more recently that complex structural and temporal patterns have been combined in a single framework [START_REF] Delvenne | Diffusion on networked systems is a question of time or structure[END_REF][START_REF] Speidel | Steady state and mean recurrence time for random walks on stochastic temporal networks[END_REF][START_REF] Hoffmann | Generalized master equations for non-Poisson dynamics on networks[END_REF]. The main purpose of this work is to go one step further along this line and to consider extreme situations where the waiting time distribution has a power-law tail with diverging average waiting times.

Diffusion with diverging average waiting times is wellknown to exhibit anomalous diffusion, with the variance * denigris.sarah@gmail.com of the process growing slower than linearly with time, when the underlying network is a regular lattice. The main novelty here is to consider a finite network with arbitrary structure, as the variance of a position is not welldefined then, and the dynamics has instead to be considered in terms of the relaxation of eigen-modes. The divergence of the average waiting time naturally leads to the formulation of the problem in terms of fractional derivatives, which already have a long tradition in the study of diffusion on regular lattices and continuous media [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF]. It is interesting to note that, despite the frequent consideration of power-law distributions and diverging moments in networks, the use of fractional derivatives is still mainly unexplored in this field of research. Rare exceptions include the recent work of Georgiou et al. [START_REF] Georgiou | Solvable non-Markovian dynamic network[END_REF] on networks with random link activation and deletion, with Mittag-Leffler interevent-time distributions and of Riascos and Mateos on the emergence of Lévy flight processes on networks through a fractional calculus approach [START_REF] Riascos | Fractional dynamics on networks: Emergence of anomalous diffusion and Lévy flights[END_REF].

This paper is organised as follows. In section 2 2.1, we introduce the model and describe the properties of its master equation. In section 2 2.2, we show that the temporal evolution can be asymptotically described by equations whose temporal derivative is fractional, which reflects the non-Markovian nature of the process. The equations reveal that the dynamical modes asymptotically relax in a power-law way, as expected, but also that the contribution of all modes remains non-negligible, for arbitrary long times. Theoretical results are tested against numerical simulations on artificial networks in section 2 2.3. Finally, we conclude and discuss the implications of our work, in particular for the problem of model reduction, in section 3.

FROM RANDOM WALKS TO FRACTIONAL DERIVATIVES

Continuous-time with power-law waiting times

We focus on the dynamics of a random walker on a network composed of N nodes. The process is encoded by structural and temporal factors. First, the N × N tran-sition matrix T determines the probability that a walker jumps from one node to one of its neighbours in the network. The transition matrix is defined as T i,j = A i,j /k j where A i,j is the adjacency matrix encoding the connections and k j is the degree,i.e. the number of connections of node j. If we consider a Markovian process, i.e. a random walk without memory, the Markov chain associated to T completely determines the statistical properties of the walker trajectories on the network.

Considering now the temporal side, we introduce the waiting time distribution ψ(t) which determines the duration of the stay on a node before the walker performs a jump. Thus this quantity directly affects the timings on the trajectory generated by T and, as a consequence some particular choices of ψ(t), the overall asymptotic behaviour shall be altered in respect to the Markovian frame. For the sake of simplicity, we have assumed in the following that the waiting time distribution is the same on each node.

Combining the two ingredients, structural and temporal, the evolution of the p i (t) ∈ R N , the probability for the walker to be at node i reads:

p p p(t) = p 0 p 0 p 0 ∞ t ψ(∆t)d∆t + t 0 ψ(t -∆t)Tp p p(∆t)d∆t, (1)
with p 0 p 0 p 0 = p p p(t = 0). In Eq. 1 the first term accounts for the possibility that a walker stays put in the node: this is clearer if we define the survival probability as

Ψ(t) = ∞ t ψ(τ )dτ = 1 - t 0 ψ(τ )dτ, (2)
which is the probability that the waiting time on a site exceeds t. The second term, on the contrary, describes the flow of probability coming from the interaction with the other nodes in the network. The master equation of this process, generalizing the classical Montroll-Weiss equation to the case of complex networks [START_REF] Montroll | Random walks on lattices II[END_REF], takes a more suitable expression in the Laplace domain

p p p(s) = 1 -ψ(s) s p 0 p 0 p 0 + ψ(s)T p p p(s) (3)
where p p p(s) is the Laplace transform of the probability p p p(s), defined as

p p p(s) = ∞ 0 p p p(t)e -st dt (4)
In the case of a Poisson process, where the waiting time is exponential by definition, this set of equations takes a simple form in the real-time domain

∂ t p p p(t) = -Lp p p(t) (5)
and one recovers the usual rate equation for continuoustime random walks, driven by the normalized Laplacian operator L = I -T , which is equivalent, up to a change of variables, to basic models for opinion dynamics [START_REF] Blondel | Convergence in multiagent coordination, consensus, and flocking[END_REF],

the position of robots in physical space [START_REF] Jadbabaie | Coordination of groups of mobile autonomous agents using nearest neighbor rules[END_REF] and synchronisation [START_REF] Strogatz | Exploring the onset of synchronization of in populations of coupled oscillators[END_REF].

In that case, the dynamics is determined by the spectral properties of the Laplacian operator. Conservation of probability implies that the dominant eigenvector, which we assume to be unique without loss of generality, has an eigenvalue zero. To determine the speed of convergence to stationarity, a key quantity is the second dominant eigenvalue, usually called the spectral gap [START_REF] Chung | Spectral Graph Theory[END_REF], which is related to important structural properties of the underlying network, such as the existence of bottlenecks and communities.

We can cast Eq. 3 in the following form to make the Laplacian operator appear:

1 ψ(s) -1 p p p(s) = 1 ψ(s) -1 1 s p 0 p 0 p 0 -Lp p p(s) (6)
The Laplacian term, from one side, entangles the dynamics carrying the interactions between nodes but, from the other side it provides the natural basis through which decompose Eq. 3 in distinct modes. Of course, since L = I -T the transition matrix and the Laplacian share the same left and right eigenvectors basis and their spectra are related by λ

(T r) α = 1 -λ (Lapl) α and λ (Lapl) α ∈ [0, 2) implies λ (T r) α ∈ (-1, 1]
. We can thus express the probability p p p(s) as a linear combination of the right eigenvec-

tors v α v α v α p p p(s) = α ãα (s)v α v α v α . (7)
Substituting this representation in Eq. 3 and projecting on the left eigenvectors basis (for the details of the calculation, see the Appendix A), we thus obtain a set of equations for the amplitudes of these modes a α

ãα (s) = 1 -ψ(s) s a init α + ψ(s)λ α ãα (s), (8)
where λ α is eigenvalue of the transition matrix associated to the

v α v α v α eigenvector Tv α v α v α = λ α v α v α v α (here and in the following λ α = λ (T r) α
, we dropped the superscript (T r) for the sake of lightness) and a init α is the projection of the initial condition p 0 p 0 p 0 = p p p(t = 0)

p 0 p 0 p 0 = α a init α v α v α v α . (9)
We obtain the formal solution which leads to exponential relaxations for Poissonian processes since the ψ(t) is exponential as well in this case. However, if the probability distribution of waiting times takes another functional form, it can entail more general relaxations, as we will display in the next Section.

Fractional dynamics

In the present work, we shall deal with waiting time distributions possessing a power law tail for t → ∞

ψ(t) ≈ β Γ(1 -β) τ β t β+1 , (10)
for 0 < β < 1 and the prefactor has been chosen to simplify further expressions. The bounds over the exponent β entail the divergence of average waiting time; nevertheless it is still possible to have an estimate of its Laplace Transform in the limit of small s by application of the Tauberian theorem [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF] 1 . Therefore its Laplace transform can be approximated by

ψ(s) ∼ 1 -(τ s) β (11)
in the limit of s → 0, which corresponds to the asymptotic regime t → ∞ we are interested in. For random walks governed by such slow decaying power tails on lattices, it is known that the average number of steps exhibits a sublinear growth, as

n(t) = 1 Γ(1 + β) t β τ β . ( 12
)
Therefore, the rate of events decreases with time, in contrast with standard processes where the average waiting is defined and n(t) ≈ t. This property tends to slow down the dynamics, leading to subdiffusion. Using Eq. 11 in the expression for the evolution of the modes in the Laplacian domain leads to

ãα (s) = (τ s) β a init α s(1 -λ α (1 -(τ s) β )) , (13)
The zero mode of the Laplacian, which is unique as we consider connected networks, corresponds to λ 1 = 1 for the transition matrix and one recovers a stationary solution, as expected:

ã1 (s) = a init 1 s (14)
For the rest of the spectrum, when λ α = 1, λ α = 0, we can cast Eq. 13 in the following form

K β α s -β + 1 ãα (s) = a init α λ α s , (15)
where K β α = 1-λα λα τ -β . In Eq. 15 we can recognize a term which carries the action of a fractional operator since we have

s -β ãα (s) = L 0 D -β t a α (t) , (16)
with D -β is the Riemann-Liouville integral operator of fractional order

β 2 0 D -β t f (t) = 1 Γ(β) t 0 f (t )dt (t -t ) 1-β . ( 17
)
We thus can interpret Eq. 15 as the Laplace transform of the integral equation

a α (t) - a init α λ α = -K β α D -β a α (t), (18)
which corresponds, applying the fractional differential operator D β on the left, to

D β a α (t) - a init α λ α = -K β α a α (t). ( 19
)
The solution to this equation is a Mittag-Leffler function [START_REF] Mathai | Special functions for applied scientists[END_REF]

a α (t) = a init α λ α E β -(k α t) β , (20)
where k α = K β α 1/β and the function is defined as

E β -(k α t) β = ∞ n=0 (-1) n (k α t) nβ Γ(βn + 1) . (21)
We observe that Eq. 15 can be interpreted as well in terms of Caputo fractional derivatives (see Appendix B). By its definition, the Mittag-Leffler function appears as a generalization of the exponential and the relaxation given by Eq. 20 approaches an exponential one in the limit α → 1. The interesting feature of these functions, in the scope of our analysis, is their asymptotic behaviour for t → ∞ which transfers the power law relaxation, proper of the waiting time distribution, to the mode themselves. Indeed, in this limit, it holds [32]

E β (-t β ) ∼ t -β Γ(1 -β) , (22)
therefore, in this limit, Eq. 20 can be approximated by

a α (t) = a init α λ α (k α t) -β Γ(1 -β) = a init α τ β 1 -λ α t -β Γ(1 -β) (23)
As for the case λ α = 0, by substitution in Eq. 13 we obtain

ãα (s) = τ β a init 0 s β-1 , (24)
which can be recast in

s -β ãα (s) = s -1 τ β a init α . (25)
With the same arguments than before, we recognize the Laplacian transform of the equation

L D -β a α (t) = L a init α τ β (26)
Thus the modes associated to the eigenvalue λ α = 0 evolve, inverting the Laplace transform, as

a α (t) = τ β a init α Γ(1 -β) t -β (27)
We can now, from Eqs. 23-27, reconstruct the asymptotic behaviour of the walker probability to be at a given node:

p p p(t) = a init 1 v 1 v 1 v 1 + N α=2 a init α τ β 1 -λ α t -β Γ(1 -β) v α v α v α . ( 28
)
To further appreciate the implications of Eq. 28, we can recast it as

p p p(t) = a init 1 v 1 v 1 v 1 + τ β Γ(1 -β)t β N α=2 a init α 1 -λ α v α v α v α , (29)
where the eigenvalues λ α in the sum are, of course, in the interval (-1, 1). It appears then evident that the spectrum of the transition matrix, and thus the Laplacian, does not organize any more the diffusion in a hierarchy of timescales, allowing to discriminate between groups of nodes where the diffusion has already relaxed to the stationary state and others which, in the worst case, shall relax with the slowest possible timescale, i.e. the second smaller eigenvalue of the Laplacian. In contrast, the introduction of a power law waiting time probability distribution, asymptotically induces in all the modes the same relaxation (Eq. 23), thus effectively wiping out the localization effect brought by the Laplacian spectrum by itself, as we shall see in Sec. 2 2.3. We note our result is coherent with the the findings of [START_REF] Delvenne | Diffusion on networked systems is a question of time or structure[END_REF], derived in situations when the waiting time distribution has an exponential tail.

Summing up the general gist of this section, we passed from the integral equation in Eq. 1 of a continuous time random walk to a fractional differential equation (Eq. [START_REF] Delvenne | Diffusion on networked systems is a question of time or structure[END_REF] in the limit of t → ∞. Now, to conclude, we would like to note that this passage might be delicate one as it is not a mere change of representation. Indeed, the choice of formulation involves the correct handling of the initial conditions: as an example, in [START_REF] Hilfer | Fractional Diffusion based on Riemann-Liouville Fractional Derivatives[END_REF], it was shown how the solution of a fractional diffusion equation can eventually turn out not being a proper probability distribution; in the present case however our starting point, Eq. 1, was shown to be correctly related to a CTRW process. Another possible caveat comes from the intrinsic non-Markovian character embedded by the fractional time derivative: as the system keeps track of its whole history, in [START_REF] Hilfer | On fractional diffusion and continuous time random walks[END_REF], the author considers two waiting time distribution, differing by a fast decaying term, that can lead to two different diffusion equations in the continuum limit; nevertheless it must be emphasized that, in the context of networks, the Laplacian operator always remains discrete when passing from a random walk description to a diffusive one and, therefore, such incongruence does not arise.

Simulations

In this Section, we turned to simulations in order to get a grip on the analytical results we have shown in the previous section. We thus explored two paths: we first considered a simulation meant to reproduce the random walk biased with a power law waiting time distribution. Through the simulation we reconstructed the evolution of the probability p p p(t) looking at the occupation of the different nodes by N walkers. The latter performed or not an hop to an adjacent node at times drawn from a waiting time distribution with various exponents β (Figs. 1ab). As displayed in Figs. 1 we indeed observe a slow contrast with the fast exponential relaxation that arises when ψ(t) is exponential as well.

We then performed a second check, which is now meant to show some implications of the coherent modes relaxation in Eq. 29: we took into account a network composed by two modules of N = 10 each that we can, according to some probability p, progressively connect (Fig. 3). Therefore for p → 1 the initial two-modules structure gets blurred and finally, there is no structural signature evidencing modularity. From the spectral point of view, this transition towards non-modularity is signalled by the spectral gap, i.e. the second smaller of the Laplacian matrix: thus, the more it approaches zero, the slower diffusion reaches the equilibrium. Therefore, in a "classical" setting in which hops are exponentially distributed there is a clear fingerprint in the dynamics entailed by a structural property. Now, in our case of study, the temporal aspect which accounts for the intermittent nature of the links, dominates over the structural one: in Fig. 4a we display the solutions for the asimptotic regime in Eq. 29 where, increasing the p parameter we progressively erase the modular structure and correspondingly the spectral gap increases, as displayed in the inset of Fig. 4a. This latter tendency of the spectral gap would give a faster relaxation to equilibrium, as in Fig. 4b, but, in our case of study, this effect is hidden by the complex temporal behaviours: in Eq. 29, we see that only two factors play a role in the relaxation, namely the initial conditions and the exponent β and, indeed, the probability distribution in Fig. 4a relaxes with the same speed on structures characterized by different spectral gaps.

CONCLUSION

In this work, we have focused on the asymptotic properties of random walk processes taking place on networks, where the distribution of waiting times between events presents a power-law tail associated to a diverging average. As shown in Sec. 2 2.1, this process can be described by a generalized master equation (Eq. 1) determining the time evolution of the probability for the walker to be in a node p p p(t). We then expressed this vector in the basis formed by the eigenvectors of the transition matrix and considered the problem in the Laplace domain, in order to derive simple algebraic equations for the relaxation of the eigenmodes of the process. As a next step, we considered the asymptotic limit t → ∞, or equivalently s → 0, where the dynamics can be rewritten in terms of Riemann-Liouville fractional operators. We first showed that all the non-stationary modes exhibit the same asymptotic power-law relaxation, independently of their eigenvalue, in contrast with the standard exponential relaxation taking place for Poisson processes. In this classical case, modes do not exhibit the same relaxation, but one observes instead a hierarchy of relaxation time-scales determined by their eigenvalue. The existence of different time scales is crucial in order to define reduced models where fast scales can be neglected in front of slower ones, but also to uncover interesting structures from the corresponding eigenvectors. In our case, however, the extreme fluctuations of the waiting times hinder the importance of the underlying structure, and imply that the amplitudes of the different eigenmodes remain of the same importance, even at large times. As expected, the power law tail of the ψ(t) also significantly slows down the relaxations, as we confirmed by means of numerical simulations.

Appendix A: Projection on the eigenmodes

The transition matrix T i,j = A ij /k j does not possess a basis of orthonormal eigenvectors due to its asymmetry. Nevertheless, its right and left eigenvector basis, {v a v a v a } and {u a } respectively, can be related with the orthonormal eigenvector basis {n a } of the symmetric matrix defined as N ij = A ij / k i k j which is related to the normalized Laplacian [START_REF] Chung | Spectral Graph Theory[END_REF]. By definition, for u a it holds i u i a T i,j = λ a u j a which can be rewritten, using the expression for T i,j , as

i u i a √ k i k j A ij k i k j = 1 k j i u i a k i N i,j = λ a u j a , (A1)
therefore the eigenvectors of the normalized laplacian can be expressed by the left eigevectors of T i,j through n j a = k j u j a . With the same reasoning one obtains the equivalent relation for the right eigenvectors n j a = v j a / k j . The usefulness of those relations appears in the attempt to separate the modes a α (t): if we consider Eq. 3 and we inject the expansions using the right eigenvectors in Eq. 7 we obtain

p i (s) = α ãα (s) k i n i = 1 -ψ(s) s α a α init k i n i α + ψ(s) α ãα (s)λ α k i n i α , (A2)
where we have expressed the right eigenvectors by the means of the normalized laplacian ones {n α }. Now, in order to isolate the different modes, we project the probability vector x on the left eigenvectors basis {u α } using, As we detailed in Sec. (B2) The differential operator acts, in the case of the Caputo definition, before the fractional integral, while the order is reversed for the R-L definition. The difference between the two thus concerns the initial conditions f 0 = f (t = 0). Indeed this fact is evident looking at their Laplace transforms which are indeed related [START_REF] Oldham | The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order[END_REF]

L C 0 D β t f (t) = L 0 D β t f (t) -s β-1 f 0 . (B3)
Therefore, coming back to Eq. 15 we can recast it as

s β ãα (s) -s β-1 a init α = λ α -1 τ β ã(s), ( B4
)
where we can recognize, on the left hand side, the Laplace transform of the Caputo derivative, obtaining

L C 0 D β t a α (t) = λ α -1 τ β L {a α (t)} , (B5)
thus illustrating, for our case of study, that the choice of the R-L derivative or the Caputo one are equivalent.

FIG. 1 .FIG. 3 .

 13 FIG. 1. Logarithmic plot of the a4(t) and a5(t) mode associated respectively to λ4 = -0.6876 and λ5 = -0.4050 for the transition matrix of the simple network in Fig. 2 for 17000 walkers and (a) β = 0.2 and (b) β = 0.4.

9 FIG. 4 .

 94 FIG. 4. (a) Relaxation towards equilibrium of the probability distribution pi at a node for different probabilities p of connection between the two modules in the case of a powerlaw waiting time distribution, here N = 10 for each module and β = 0.7. (inset) Variation of the spectral gap for the Laplacian matrix spectrum λ L = 1 -λ T in correspondence to the bridging probability p. (b) Probability distribution p1(t)when ψ(t) ∝ exp(-t/τ ) for the same networks of (a). Different eigenvalues now correspond to different relaxation times, and modes associated to small eigenvalues exhibit a slower relaxation, therefore dominating for long times.

  2 2.2, Eq. 15 admits an interpretation as the Laplace Transform of a fractional differential equation in terms of Riemann-Liouville (R-L) operators. Equivalently, it might be interpreted using the Caputo definition of such operators. The Caputo fractional derivative C 0 D β t reads (in what follows we restrict ourselves to the case 0 < β < 1)C 0 D β t (f (t)) = I 1-β (f (t)dt (tt ) -β ,(B1) while the R-L fractional derivative 0 D β t is defined as0 D β t (f (t)) = dt (tt ) -β .

  again, their relation with the {n α }

	i	n i β √ k i	p i (s) =	α	ãα (s)	i	n i β √ k i	k i n i α = a β (s) =
			1 -ψ(s) s	α	a α init	i	n i β √ k i	k i n i α +
			ψ(s)	α	ãα (s)λ α	i	n i β √ k i	k i n i α . (A3)
	Since i n i β n i α = δ α,β the previous equation gives
		a β (s) =	1 -ψ(s) s	a β init + ψ(s)λ β a β (s),	(A4)
	which is our Eq. 8.

Appendix B: Connection with Caputo fractional operators

More precisely, we applied the Hardy-Littlewood Tauberian theorem to Ψ(t) in Eq

. 2 in the form presented in[START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF], Ch. XIII Sec. 5, p.446, Theorem 4.

In what follows, the bounds of integration shall always be 0 and t, thus we shall omit them to have a lighter notation.