## Abstract

We present a class of continuous-time random walks (CTRWs), in which random jumps are separated by random waiting times. The novel feature of these CTRWs is that the jumps are clustered. This introduces a coupled effect, with longer waiting times separating larger jump clusters. We show that the CTRW scaling limits are time-changed processes. Their densities solve two different fractional diffusion equations, depending on whether the waiting time is coupled to the preceding jump, or the following one. These fractional diffusion equations can be used to model all types of experimentally observed two power-law relaxation patterns. The parameters of the scaling limit process determine the power-law exponents and loss peak frequencies.

## 1. Introduction

The aim of this paper is to show that a compound subordination approach to anomalous diffusion, based on clustered continuous-time random walk (CTRW) methodology, provides useful tools to study relaxation phenomena in complex systems. This approach allows us to identify stochastic origins of all parameters describing the non-Debye fractional power-law relaxation behaviour observed in complex materials, such as supercooled liquids, amorphous semiconductors and insulators, polymers, disordered crystals, molecular solid solutions, glasses and so on (Jonscher 1983; Scher *et al.* 1991; Havriliak & Havriliak 1994). Relaxation data, collected in both time and frequency domains, show that the dielectric susceptibility *χ*(*ω*)=*χ*^{′}(*ω*)−*iχ*^{′′}(*ω*) exhibits the following asymptotic dependence on frequency:
1.1
where Δ*χ*(*ω*)=*χ*(0)−*χ*(*ω*), the exponents *n* and *m* fall in the range (0,1), and *ω*_{p} denotes the loss peak frequency, the reciprocal of a time constant characteristic for a given material. Experimental evidence (Havriliak & Havriliak 1994; Jonscher 1996) shows that relaxation data fall into two major groups: typical relaxation behaviour, for which *m*≥1−*n*, and less typical relaxation, where *m*<1−*n*. The boundary case, *m*=1−*n*, is known as Cole–Cole (CC) relaxation.

To fit measured dielectric susceptibility data, the empirical Havriliak–Negami (HN) function
1.2
has been proposed (Jonscher 1983; Havriliak & Havriliak 1994; Jonscher 1996). For *γ*=1, the formula (1.2) describes CC relaxation. The original HN function satisfies the power-law properties (1.1) with *n*=1−*αγ* and *m*=*α*, characteristic of typical relaxation behaviour. To model less typical power-law dielectric data, an extended range of the power-law exponents, 0<*α*<1 and 0<*αγ*<1, has been proposed for the HN function (Havriliak & Havriliak 1994). However, the derivation of the HN function (1.2), based on the fractional Fokker–Planck equation (Kalmykov *et al.* 2004) and the CTRW (Weron *et al.* 2005; Jurlewicz *et al.* 2008), does not lead to values of *γ*>1. Only recently, by employing the subordination approach (Stanislavsky *et al.* 2010; Weron *et al.* 2010), developed in the last decade in the theory of anomalous diffusion (Sokolov 2001, 2002; Meerschaert *et al.* 2002; Meerschaert & Scheffler 2004; Piryatinska *et al.* 2005; Magdziarz & Weron 2006; Sokolov & Klafter 2006; Magdziarz *et al.* 2007, 2008; Lubelski *et al.* 2008), has the effective picture underlying all typical and less typical relaxation patterns been found. In contrast to earlier studies (Weron & Kotulski 1996; Metzler *et al.* 1999) where the anomalous-diffusion approach (based on a decoupled CTRW) has led to CC relaxation only, the compound subordination framework allows one to derive rigorously (Stanislavsky *et al.* 2010; Weron *et al.* 2010) the original HN function as well as the novel version
1.3
termed Jurlewicz, Weron and Stanislavsky (JWS) by Trzmiel *et al.* (2011). The JWS function (1.3) satisfies the power-law properties (1.1) with *n*=1−*α* and *m*=*αγ* fulfilling the relation *m*<1−*n*, and hence it is appropriate for description of less typical relaxation data.

In this paper, we start by describing the clustered CTRW, accounting for interactions between space jumps and waiting times. Then, we examine the relaxation behaviour (1.1) with any *n* and *m* falling into the range (0,1). We derive the governing diffusion equations and the time-domain relaxation functions corresponding to the HN and JWS shape functions, equations (1.2) and (1.3), respectively. In our study, we draw special attention to the origins of the loss peak frequency *ω*_{p}. The analysis moves forward to a linkage between kinetics and thermodynamics in relaxing materials. The point is that the theory of CTRW processes has free parameters (diffusive coefficients, indices of random processes and others) which, as we show, can depend on temperature. This feature is well known in experimental physics, for example, as the Arrhenius law. All this opens a way to relate different branches of physics.

## 2. Continuous-time random walk and anomalous diffusion

The CTRW process *W*(*t*) determines the location reached at time *t*, for a particle performing a random walk, in which random particle jumps are separated by random waiting times (Klafter & Sokolov 2011). Together, the jumps and the inter-jump waiting times form a sequence of independent and identically distributed random vectors (*R*_{i},*T*_{i}), *i*≥1. Both the length and the direction of the random jump variable *R*_{i} can depend on the waiting time *T*_{i}. The position of the particle after *n* jumps reads
2.1
In the classical waiting-jump CTRW idea of Montroll & Weiss (1965), in which the jump *R*_{i} occurs after the waiting-time, the random number of particle jumps performed by time *t*>0 is given by the renewal counting process
2.2
where
2.3
is the time of the *n*th jump. Then, the location of a particle at time *t* is given by the random sum
2.4
In the alternative jump waiting CTRW scenario, the particle jump *R*_{i} precedes the waiting time. Now, the counting process *N*_{t}+1 gives the number of jumps by time *t*, and the particle location at time *t* is given by
2.5
This is called the over-shoot CTRW, or briefly OCTRW (Jurlewicz *et al.* in press). In summary, the CTRW process *W*^{−}(*t*) and the OCTRW process *W*^{+}(*t*) are obtained by the subordination of the random walk *R*(*n*) to the renewal counting process *N*_{t}, and the first passage process *N*_{t}+1, respectively. Let us mention that, in general, subordination modifies a random process, replacing the deterministic time index by a random clock process, which usually represents a second source of uncertainty (see Feller (1971) or Sato (1999) for more details).

As far as the diffusion mechanism underlying the relaxation phenomenon is concerned, the diffusion limit (diffusion front) of the random sum given by (2.4) or (2.5) should be considered. The diffusion front of a process *X*(*t*) is represented by the asymptotic behaviour of the rescaled process *f*(*c*)*X*(*ct*) when the dimensionless time-scale coefficient *c* tends to infinity and the space-rescaling function *f*(*c*) is chosen appropriately. The possible general form of the diffusion limit, as well as conditions providing its existence, have already been precisely determined for both CTRW and OCTRW processes (Weron *et al.* 2010; Jurlewicz *et al.* in press). When the jumps *R*_{i} and the waiting times *T*_{i} are stochastically independent (i.e. uncoupled), the CTRW and OCTRW diffusion limits and, as a consequence, the corresponding type of relaxation, are the same. On the other hand, if the coupled case (i.e. dependent coordinates in the random vector (*T*_{i},*R*_{i})) is considered, the waiting-jump and jump-waiting schemes may lead to essentially different relaxation patterns.

An important and useful example of coupled CTRW, different from the most popular Lévy walk (Klafter & Sokolov 2011), has been identified using the clustered CTRW concept, introduced by Weron *et al.* (2005) and developed further by Jurlewicz *et al.* (2011). While in the Lévy walk, the jump size is fully determined by the waiting time (or equivalently, by flight duration), in a clustered CTRW, coupling arises from clustering of a random number of jumps (figure 1). The clustered CTRW scheme is relevant to numerous physical situations, including the energy release of individual earthquakes in geophysics, the accumulated claims in insurance risk theory, and the random water inputs flowing into a reservoir in hydrology (Weissmann *et al.* 1989; Klafter & Zumofen 1994; Huillet 2000). In all these cases, summing the individual contributions yields the total amount (in general, random) of the studied physical magnitude over certain time intervals. In the clustered CTRW, both the waiting time and the subsequent jump are random sums with the same random number of summands, and as a consequence, this type of CTRW is coupled, even if the original CTRW before clustering had no dependence between the corresponding waiting times and jumps. If the random number of jumps in a cluster has a heavy-tailed distribution, then the effect of clustering on the limiting distribution can be profound. In this case, the OCTRW jump-waiting scheme and the traditional CTRW waiting-jump model are significantly different in both their diffusion limit, and their governing equations (Weron *et al.* 2005; Jurlewicz *et al.* 2011).

Considering those cases crucial for modelling of relaxation phenomena, assume that the waiting times *T*_{i} have a heavy-tailed distribution with parameter 0<*α*<1; i.e. for some *τ*_{0}>0, we have for large *t*. Moreover, let the jump distribution be symmetric and belong to the normal domain of attraction of a symmetric Lévy-stable law with the index of stability *η*, 0<*η*≤2; i.e. for some *ρ*_{0}>0, we have for large *x* if 0<*η*<2 or , if *η*=2. Then, a heavy-tailed distribution of cluster sizes with tail exponent 0<*γ*<1 yields different anomalous diffusion limits
2.6
2.7
of clustered CTRW and clustered OCTRW, respectively (Weron *et al*. 2010; Jurlewicz *et al.* 2011). Here, *C* is a positive constant dependent on the tail exponents *α*, *η* and proportional to the scaling parameter *ρ*_{0}. Each limit takes the form of the parent process , coming from the diffusion limit of the cumulative-jump sequence (2.1), subordinated to a pair of time changes. The inner time change is an inverse-time *α*-stable subordinator that captures the effect of the waiting-time sequence. The novelty of the clustered CTRW limit is the intermediate time change, , , that reflects the clustering, for under- and over-shooting processes.

The processes , and in the representation (2.6) of the clustered CTRW diffusion front are stochastically independent (as well as those in (2.7)), since the clustered CTRW (and OCTRW) models assume stochastic independence between the waiting times, jumps and cluster sizes. The parent process is a symmetric *η*-stable Lévy process (in particular, for *η*=2, it is just a standard Brownian motion). The directing process is defined as
where is a strictly increasing stable Lévy process with the stability index *α* being the diffusion front of the cumulative-waiting-time sequence (2.3). The stable Lévy process with the stability index *η* is 1/*η*-self-similar, and hence (Embrechts & Maejima 2002). Similarly, is 1/*α*-self-similar, which implies the inverse scaling for the inverse subordinator (Meerschaert & Scheffler 2004). To ease notation, and without any real loss of generality, we assume henceforth that and .

The under-shooting and over-shooting processes have the same distributions as the processes and , respectively, for the value of parameter *α* equal to *γ*, where is the left-limit process corresponding to . The names refer to the property that, for any time *t*≥0,
2.8
so that the under- and over-shooting processes under- and over-estimate the exact instant of time *t*≥0, respectively. As shown in Weron *et al.* (2010) and Jurlewicz *et al.* (2011), both processes are 1-self-similar and, for a fixed *t*≥0, the distribution of is the same as *tJ* for the random variable *J* with the probability density function (PDF)
2.9
while , where the PDF of *Z* reads
2.10
The functions *p*_{J}(*x*) and *p*_{Z}(*x*) are easily recognized as special cases of the well-known beta densities. The PDF *p*_{J}(*x*) concentrates near 0 and 1, which means that the most probable values for occur near 0 and *t*. Moreover, has finite moments of any order that can be calculated directly from the density (2.9). On the other hand, *p*_{Z}(*x*) concentrates near 1 and possesses the power-law property *p*_{Z}(*x*)∝*x*^{−γ−1} for large *x*. Hence, the values of not only concentrate near *t*, but also are dispersed along the positive half-line, and even the first moment of diverges because of the heavy tail of its distribution.

Assuming that the simple normalizing function *f*(*c*)=*c*^{−α/η} was used to arrive at the diffusion limits, the positive constant parameter *C* in (2.6) and (2.7) can be explicitly computed,
2.11
where
Thus, *C* is just the scaling jump parameter *ρ*_{0} multiplied by a proportionality constant that only depends on the tail exponents *α* and *η*. In particular, the constant *C* is not influenced by the clustering parameters.

## 3. Two power-law responses and anomalous diffusion

From rather general assumptions, the theoretical attempt to model non-exponential relaxation phenomena is based on the idea of relaxation of an excitation undergoing diffusion in the system under study (Metzler & Klafter 2000). Consequently, the time-domain relaxation function *ϕ*(*t*) can be related to the dielectric susceptibility *χ*(*ω*) as
In the framework of a one-dimensional CTRW, the time-domain relaxation function *ϕ*(*t*) is given in terms of the temporal decay of a given mode *k* by the Fourier transform
3.1
where *k*>0 has the physical sense of a wavenumber, and denotes the diffusion front under consideration. If experimental tools probe the system for a given mode, the above formula gives the temporal relaxation response after a macroscopic excitation. The relaxation patterns *ϕ*(*t*) are determined by the stochastic properties of the jumps, the inter-jump times and the detailed construction of the counting process.

For the cluster CTRW and OCTRW models, the anomalous diffusion fronts in (2.6) and (2.7) have frequency-domain dielectric susceptibility functions that can be identified as the JWS and HN functions, equations (1.3) and (1.2), respectively, see Stanislavsky *et al.* (2010) and Weron *et al.* (2010). The characteristic material constant *ω*_{p}, appearing in both functions, takes the form
3.2
where *C* is defined by (2.11), reflecting its stochastic origins. Representations (2.6) and (2.7) allow us to develop analytical formulae for the corresponding time-domain relaxation functions *ϕ*_{JWS}(*t*) and *ϕ*_{HN}(*t*) in terms of the three-parameter Mittag–Leffler function (Mathai *et al.* 2009),
Here, is Appell’s symbol with (*γ*,0)=1. Observe that
3.3
where *p*^{±}(*x*,*t*) are PDFs of the anomalous diffusion processes . Hence, the Fourier–Laplace (FL) images of the functions *p*^{±}(*x*,*t*) are just the Laplace images of the corresponding time-domain JWS and HN relaxation functions. According to the results of Jurlewicz *et al.* (2011) and Stanislavsky & Weron (2010), we can derive *p*^{−}(*x*,*t*) as a mild solution^{1} of the following fractional pseudo-differential equation:
where *δ*(*x*) is the Dirac delta function. Equivalently,
3.4
and hence
3.5
with *ω*_{p} as in (3.2). Using the relation
we obtain
3.6

Similarly, we obtain that *p*^{+}(*x*,*t*) is a mild solution of
where *p*_{R}(*x*,*t*) and *p*_{S}(*τ*,*t*) are the PDFs of and , respectively (Stanislavsky & Weron 2010; Jurlewicz *et al.* 2011). This implies that
3.7
As a consequence, we have
3.8
3.9
Plots of the HN and JWS time-domain relaxation functions, obtained by numerical simulations, are presented in figure 2. The short- and long-time behaviours of these functions, which exactly follow the high- and low-frequency power laws (1.1), read
Figures 3 and 4 illustrate power-law properties of the HN and JWS time-domain relaxation functions.

## 4. Temperature dependence of free parameters

In the theory of relaxation (Scher *et al.* 1991; Jonscher 1996; Magdziarz & Weron 2006), considerable effort has gone into finding any function that could scale such quantities as relaxation time, viscosity and diffusion coefficient. The scaling idea provides a linkage between thermodynamics and kinetics in relaxing materials near the phase (glass-like) transition. For the clustered CTRW model of anomalous diffusion presented here, the parameters *ρ*_{0}, *τ*_{0}, *α*, *γ* and *η* are free, in the sense that they are independent of time and space coordinates. On the other hand, they may vary with the changing thermodynamic state of the physical system (characterized by temperature or pressure). Their temperature/pressure dependence is given by transition-state theory.

The indices *α* and *γ* directly determine the power-law exponents *m* and *n* characterizing the asymptotic behaviour of the relaxation patterns in time and in frequency. The parameter *α* is associated with the temporal properties of the relaxing system under consideration. The inverse *α*-stable process accounts for the resting time of a walker (Baeumer *et al.* 2005), and the walker randomly moves all the time only if *α* = 1, whereas it rests forever when *α*=0. The index *γ* describes the spatial-temporal clustering present in the relaxing system. If *γ*=1, the influence of clustering on the diffusion limit disappears, whereas the case 0<*γ*<1 corresponds to a more complicated structure, including the over- or under-shooting component of the compound time change. If *γ* is close to 0, large clusters of jumps separated by long inter-jump times become more and more dominant. In contrast, as *γ* tends to 1, the clustered spatial-temporal steps become comparable to the steps without clustering.

Changes of temperature and pressure can modify the coefficients *α* and *γ*, as in the case of so-called beta relaxation (Jonscher 1983, 1996), as well as the tail exponent *η* of the random jumps, which characterizes the diffusion dynamics. It is surprising that *α* and *γ* are independent of the thermodynamic quantities, in the case of so-called alpha relaxation (Jonscher 1983, 1996). If *η* is also temperature/pressure independent, we can explain the temperature/pressure dependence of the loss peak frequency
(see equations (3.2) and (2.11)) in terms of the scaling coefficients *ρ*_{0} and *τ*_{0}. The spatial scaling coefficient *ρ*_{0} should increase with temperature, making long jumps more likely (since for large jump length *x*, so it increases with *ρ*_{0}). In contrast, the temporal scaling coefficient *τ*_{0} should decrease, resulting in shorter waiting times. These properties follow the classical Arrhenius law, assuming an exponential character of *ρ*_{0} and 1/*τ*_{0} depending on the pressure and the reciprocal of temperature. The loss peak frequency *ω*_{p} is thus proportional to , where *k*_{B} is the Boltzmann constant, *T* the temperature, *P* the pressure and *E* and *V* are the activation energy and volume, respectively.

The proposed stochastic model for two power-law relaxation can also reveal the interrelation between the empirical thermodynamic behaviour of the relaxing system and the anomalous diffusion parameters, when *α*, *γ* and *η* are also dependent on temperature/pressure. However, this analysis requires more advanced techniques, and is beyond the scope of this paper.

## 5. Conclusions

In this paper, we have shown that clustered CTRWs and their diffusion limits illuminate the role of random processes in the parametrization of relaxation phenomena. The key result is that the empirical relaxation functions can be derived from diffusion models, introducing free scaling parameters, independent of time and space coordinates. This allows one to consider the parameters as material coefficients, depending on thermodynamic quantities. As we have demonstrated, the *η*-stable parent process that models particle jumps does not change the type of the relaxation response, and it only affects the material constant *ω*_{p} by determining the spatial features of the anomalous diffusion *W*^{±}(*t*). The index *α* of the process that codes the waiting times between jumps, and the index *γ* of the clustering process determine the power-law behaviour of the relaxation function in time and frequency. These coefficients characterize the complex dynamics of relaxing systems.

## Acknowledgements

The work of K.W. and A.J. was partially supported by the Ministry of Sciences and Higher Education project PB NN 507503539, The work of M.M.M. was partially supported by USA National Science Foundation grants DMS-1025486 and DMS-0803360, and USA National Institutes of Health grant R01-EB012079-01. A.S. is grateful to the Institute of Physics and the Hugo Steinhaus Center for Stochastic Methods for pleasant hospitality during his visit to Wrocław University of Technology.

## Footnotes

↵1 A mild solution to a differential equation is a function whose Laplace or Fourier transform solves the corresponding algebraic equation. See Baeumer

*et al.*(2005) or Meerschaert & Scheffler (2008) for more details.

- Received November 28, 2011.
- Accepted January 5, 2012.

- This journal is © 2012 The Royal Society