Stochastic Tools Hidden Behind the Empirical Dielectric Relaxation Laws

TL;DR

This paper reviews stochastic modeling techniques for dielectric relaxation in disordered systems, highlighting their ability to explain universal kinetic patterns across various complex materials.

Contribution

It introduces recent stochastic tools and probabilistic models that analyze the relaxation phenomena and their universal features in complex disordered systems.

Findings

Stochastic models explain universal dielectric relaxation patterns.

Diffusion and tempered processes provide insights into anomalous dynamics.

Relationships among stochastic approaches enhance understanding of relaxation mechanisms.

Abstract

The paper is devoted to recent advances in stochastic modeling of anomalous kinetic processes observed in dielectric materials which are prominent examples of disordered (complex) systems. Theoretical studies of dynamical properties of those "structures with variations" [Goldenfield and Kadanoff 1999 Science 284 87--9] require application of such mathematical tools by means of which their random nature can be analyzed and, independently of the details differing various systems (dipolar materials, glasses, semiconductors, liquid crystals, polymers, etc.), the empirical universal kinetic patterns can be derived. We begin with a brief survey of the historical background of the dielectric relaxation study. After a short outline of the theoretical ideas providing the random tools applicable to modeling of relaxation phenomena, we present probabilistic implications for the study of the relaxation-rate distribution models. In the framework of the probability distribution of relaxation rates we consider description of complex systems, in which relaxing entities form random clusters interacting with each other and single entities. Then we focus on stochastic mechanisms of the relaxation phenomenon. We discuss the diffusion approach and its usefulness for understanding of anomalous dynamics of the relaxing systems. We also discuss extensions of the diffusive approach to the systems under tempered random processes. Useful relationships among different stochastic approaches to the anomalous dynamics of complex systems allow us to get a fresh look on this subject. The paper closes with a final discussion on achievements of stochastic tools describing the anomalous time evolution of the complex systems.

Tables

Table 1. Table 1: Relationship between internal random properties of clustered complex systems and their relaxation responses. The indices of the column “Power-law exponents” correspond to the values 0 < α , γ ≤ 1 formulae-sequence 0 𝛼 𝛾 1 0<\alpha,\gamma\leq 1 in Eqs.( 63 ) and ( 66 ), m 𝑚 m and n 𝑛 n are related to Eqs.( 3 ) and ( 4 ).

	Power-law
Law	exponents		$N_i$	$M_j$	${\beta }_{iN}$
	$m$	$1-n$
D	1	1
				long tail
CD	1	$\gamma$		$\gamma$
			long tail		long tail
CC	$\alpha$	$\alpha$	$\alpha$		$\alpha$
mirror				long tail
CD	$\gamma$	1		$\gamma$
			long tail	long tail	long tail
HN	$\alpha$	$\alpha \gamma$	$\alpha$	$\gamma$	$\alpha$
			long tail	long tail	long tail
JWS	$\alpha \gamma$	$\alpha$	$\alpha$	$\gamma$	$\alpha$

Table 2. Table 2: Relationship between the internal random properties of clustered complex systems and their fractional short-time power-law relaxation responses. The indices of the column “Short-time power-law exponent” corresponds to Eqs. ( 7 ), and 0 < α , γ ≤ 1 formulae-sequence 0 𝛼 𝛾 1 0<\alpha,\gamma\leq 1 are related to Eq.( 73 ).

	Short-time
Law	power-law	$N_i$	$M_j$	${\beta }_{iN}$
	exponent $n$
			long tail
CD	$1-\gamma$		$\gamma$
KWW	$1-\alpha$			long tail
				$\alpha$
			long tail	long tail
GG	$1-\alpha \gamma$		$\gamma$	$\alpha$

Table 3. Table 3: Relation between the parameters α 𝛼 \alpha , γ 𝛾 \gamma and σ 𝜎 \sigma in Eq.( 149 ) leading to particular forms of relaxation fitting functions.

		$\sigma \ge 0$	Type	Ref.
any	any	any	given by (149)	[83]
1	1	any	D	[20]
any	1	0	CC	[23, 24]
any	1	1	CD	[25, 26]
any	1	any	pseudo CD	[142]
any	any	0	JWS	[80]
any	$\alpha$	any	pseudo CC	–
1	any	any	mirror CD	–

Table 4. Table 4: Survival probability of the initial non-equilibrium state of a system in different probability approaches.

Probability	Survival probability of	Key notations
approach	the non-equilibrium state
master
equation	$\varphi \left(t\right)=e^{-{\int }_0^{t}r\left(s\right)𝑑s}$	Eq.(18)
relaxation
rate	$\varphi \left(t\right)=⟨e^{-\tilde{\beta }t}⟩$	Eq.(31)
distribution
two-state	$\varphi \left(t\right)=⟨e^{-{\omega }_p\tau }⟩$	Eq.(90)
relaxation
diffusion	$\varphi \left(t\right)=\{\begin{array}{cc}⟨e^{ik\tilde{R}\left(t\right)}⟩& \\ ⟨e^{-k\overline{R}\left(t\right)}⟩& \end{array}$	Eqs.(95), (96)

Equations

$${\varphi }_{\mathrm{HN}}^{\ast }(\omega )=\frac{1}{[1+(i\omega /{\omega }_p{)}^{\alpha }{]}^{\gamma }},0<\alpha ,\gamma \le 1,$$

$${\varphi }_{\mathrm{KWW}}(t)=e^{-(t/{\tau }_p{)}^{\alpha }}$$

$$\chi (\omega )\propto {\left(\frac{i\omega }{{\omega }_p}\right)}^{n-1}\mathrm{for}\omega \gg {\omega }_p,$$

$$\Delta \chi (\omega )\propto {\left(\frac{i\omega }{{\omega }_p}\right)}^m\mathrm{for}\omega \ll {\omega }_p,$$

$${\varphi }^{\ast }(\omega )={\int }_0^{\infty }{e}^{-it\omega }d(1-\varphi (t))$$

$${\varphi }^{\ast }(\omega )={\int }_0^{\infty }{e}^{-it\omega }f(t)dt$$

$$f(t)\propto \{\begin{array}{ll}(t/{\tau }_p{)}^{-n}& \mathrm{for}t\ll {\tau }_p\\ V(t)& \mathrm{for}t\gg {\tau }_p\end{array}$$

$$V(t)=\{\begin{array}{ll}(t/{\tau }_p{)}^{-m-1}& \mathrm{for}\mathrm{HN}\mathrm{and}\mathrm{CC}\mathrm{relaxation},\\ \exp \left[-(t/{\tau }_p{)}^{\alpha }\right]& \mathrm{for}\mathrm{KWW}\mathrm{relaxation},\\ \exp \left(-t/{\tau }_p\right)& \mathrm{for}\mathrm{CD}\mathrm{relaxation}.\end{array}$$

$$F(x)={\int }_{-\infty }^{x}f(y)dy$$

$$⟨X^n⟩={\int }_{-\infty }^{\infty }{x}^{n}f(x)dx.$$

$$⟨g(X)⟩={\int }_{-\infty }^{\infty }g(x)f(x)dx.$$

$$F(a_{1}x+b_1)\ast F(a_{2}x+b_2)=F(ax+b)$$

$$F_1\ast F_2=\int F_1(x-y)d{F}_2(y).$$

$$a=(a_1^{\alpha }+a_2^{\alpha }{)}^{1/\alpha }\mathrm{with}0<\alpha \le 2$$

$$h(s)={\int }_{-\infty }^{\infty }\exp (isx)dF(x).$$

$$\log h(s)=\{\begin{array}{l}i\gamma s-\sigma |s{|}^{\alpha }\{1-i\beta \mathrm{sgn}(s)\tan (\pi \alpha /2)\}\text{mbox}if\alpha \ne 1,\\ i\gamma s-\sigma |s|-i\beta (2/\pi )\sigma s\log |s|\text{mbox}if\alpha =1,\end{array}$$

$$\underset{x\to \infty }{\lim }\frac{\mathrm{Pr}(X>x)}{\sigma x^{-a}}=1$$

$$\int f_X(x/y)f_Y(y)\frac{dy}{y}.$$

$$\int f_X(xy)f_Y(y)ydy.$$

$$\int f_X(x-y)f_Y(y)dy.$$

$${\varphi }_{\mathrm{D}}(t)=\exp (-t/{\tau }_p)$$

$$\varphi (t)=\sum _{i=1}^n{p}_i\exp \left(-t/{\tau }_i\right),$$

$$\sum _{i=1}^n{p}_i=1.$$

$$\varphi (t)={\int }_0^{\infty }\exp (-t/\tau )w(\tau )d\tau ,$$

$$\varphi (t)=⟨e^{-t/\tilde{T}}⟩$$

$$\frac{d\varphi (t)}{dt}=-r(t)\varphi (t),\varphi (0)=1.$$

$$p(t,dt)=\Pr\bigl{(}t\leq\tilde{\theta}\leq t+dt|\tilde{\theta}\geq t\bigr{)}\,,$$

$$p(t,dt)=-\frac{\mathrm{Pr}(\tilde{\theta }\ge t+dt)-\mathrm{Pr}(\tilde{\theta }\ge t)}{\mathrm{Pr}(\tilde{\theta }\ge t)},$$

$$p(t,dt)=-d\ln\Pr\bigl{(}\tilde{\theta}\geq t\bigr{)}\,.$$

$$r(t)=-\frac{d}{dt}\ln \mathrm{Pr}(\tilde{\theta }\ge t)$$

Full text

1. 1 Introduction

2. 1.1 Origins of the theory of relaxation

3. 1.2 Experimental peculiarity of relaxation

4. 2 Definitions and terminology

5. 2.1 Limit theorems

6. 2.2 Lévy $\alpha$-stable distributions

7. 2.3 Mixtures of distributions

8. 2.4 Stochastic processes. Subordination

9. 3 Relaxation function as the initial-state survival probability. Probabilistic models

10. 3.1 Transition probability

11. 3.2 Microscopic stochastic scenario

12. 3.3 Stretched exponential relaxation

13. 3.4 Conditionally exponential decay model

14. 3.5 Relaxation of hierarchically clustered systems

15. 3.5.1 Havriliak-Negami function

16. 3.5.2 Infinite mean cluster sizes

17. 3.5.3 Finite mean cluster sizes

18. 3.6 Spatial randomness

19. 3.7 Probabilistic vs deterministic modeling

20. 4 Relaxation in two-state systems

21. 4.1 Stochastic arrow of time

22. 4.2 Memory function

23. 4.3 Memory function formalism for empirical laws of relaxation

24. 5 Anomalous diffusion approach. Stochastic processes

25. 5.1 Uncoupled continuous time random walk. Diffusion. Cole-Cole relaxation

26. 5.2 Coupling between the very large jumps in physical and operational times

27. 5.3 Anomalous diffusion with under- and overshooting subordination

28. 5.4 Relaxation process from compound subordinators

29. 5.5 Clustered continuous time random walk. Compound counting process

30. 6 Extensions

31. 6.1 Simple tempered relaxation

32. 6.2 Tempered relaxation with clustering patterns

33. 7 Subordination impact in understanding of relaxation scenarios

34. 8 Outline

Dedicated in memory of Leo Kadanoff

1 Introduction

Motion of charges, their accumulation and discharge are the basis of many physical, chemical and biological processes in nature. Undoubtedly, many-body interactions [1, 2, 3] play an appreciable role in the time evolution of such systems. Besides, the systems themselves are weakly or strongly disordered. This aspect is very versatile. Defects, vacancies and dislocations are frequently present in real materials [4]. Amorphous materials possess a marked departure from crystalline order [5], and a perfect (ideally ordered) crystal is difficult to find in nature. There exists a great variety of materials that have a local order in few atoms or molecules, but their structure becomes disordered on larger length scales. In a consequence, these effects induce relaxation processes inseparably linked with disorder in the systems. The relaxing entities – dipoles, traps, ions and so on, interact not only among themselves, but also with the surrounding medium to modify disorder of this medium and to affect other entities. The transformations include time fluctuations in potentials seen by each entity and essentially act as a noise source. On the other hand, they form a complex potential landscape with many local minima separated by barriers of all scales, trapping and untrapping the entity orbits in a self-similar hierarchy of cantori. Consequently, the motion of entities can be very similar to a random walk. It is not surprising that parallels, suggested in literature [6, 7, 8, 9], to be drawn between relaxation and diffusion.

The relaxation properties of various complex systems (amorphous semiconductors and insulators, polymers, molecular solid solutions, glasses, etc) have attracted an immediate interest of scientists and technologists for a long time [10] (and the references therein). This gave a huge wealth of experimental data. The data analysis discovered the “universality” of relaxation patterns [11] per se that is enclosed into fractional power laws of relaxation responses (in frequency and time) for a very wide range of materials. The fascinating behavior covers 17 decades ($\sim$10*-5*-1012 Hz in frequency or $\sim$10*-12*-105 s in time), and a theoretical explanation of the universal relaxation response is one of the most difficult problems of Physics today, as any realistic physical treatment of relaxation has to take into account the stochastic or probabilistic representation of the system’s behavior. Most of the interpretations in literature (see their comprehensive discussion, for example, in [10, 11]) explain only a limited number of characteristics of the relaxation processes in complex systems. Without any doubt the physical processes going on in disordered media are complex and independent on the details of the systems under investigation. Therefore, any simple interpretation based on one or two observational facts will not explain all the features of relaxation patterns self consistently. This review is just devoted to recent advances in the theory of relaxing complex media. Our approach is based on limit theorems of probability theory. The usefulness of the theorems is that they allow us to connect microscopic stochastic dynamics of relaxing entities with the macroscopic deterministic behavior of the systems as a whole. The universal macroscopic relaxation response appears not to be attributed to any particular entity taken from those forming such a system. Any excited complex system tending to equilibrium passes from less disordered states to more disordered ones. Macroscopic evolution of the system is a result of averaging over local random properties of system’s entities. The complex system is a marvelous “supercomputer” capable for the procedure, and the limit theorems of probability theory play the role of its “software”.

The complex systems and the investigation of their structural and dynamical properties have established on the physics agenda almost three decades ago. These “structure with variations” [1] are characterized through (i) a large diversity of elementary units, (ii) strong interactions between the units, (iii) a non-predictable or anomalous evolution in course in time [12]. Their study play a dominant role in exact and life sciences, including a richness of systems such as glasses, liquid crystals, polymers, proteins, biopolymers, organisms or even ecosystems.

1.1 Origins of the theory of relaxation

Experimental and theoretical studies of relaxation phenomena have a long history. The first measurements of electrical relaxation were carried out for alkali ions in the Leyden jar (a glass) in 1847 and 1854 by R. Kohlrausch [13, 14], and the observations of mechanical relaxation in the natural polymer, silk, in 1863 and 1866 were continued by his son, F. Kohlrausch [15]. The concept of “relaxation time” into physics and engineering was introduced by J. C. Maxwell in 1867 [16]. As was shown by M. J. Curie [17] and E. von Schweidler [18] the dielectric relaxation response in the time domain can be described by a short-time power-law dependence. Perhaps, P. Debye in 1913 was the first who derived soundly the relaxation relationship based on principles of statistical mechanics [19, 20]. For this purpose he used Einstein’s theory of Brownian motion [21, 22] to consider the collisions between a rotating dipolar molecule and its neighboring other molecules in the liquid under the assumption that the only electric field acting on the molecule is an external field. Consequently, the Debye relaxation law is expressed in terms of rotational Brownian motion, and it has an exponentially decaying form in time domain. The physical mechanism underlying the Debye law is obviously simpler than the one underlying the stretched exponential relaxation found by Kohlrausch [14]. Rapid developments of science and technology in the twentieth and twenty-first centuries created a wide variety of new materials with non-exponential behavior in relaxation properties. Over the past 100 years, many empirical relaxation laws, which are regarded as generalizations of the Debye (D) relaxation law, have been discovered. Among the most known, and frequently utilized to analyze the frequency domain measurements, are the Cole-Cole (CC) law (1941-1942) [23, 24], the Cole-Davidson (CD) law (1950-1951) [25, 26] and the Havriliak-Negami (HN) law (1966-1967) [27, 28]. In fact, all they were established empirically. Their form is convenient to write as

[TABLE]

$\alpha=\gamma=1$ being the D case; $\alpha=1,0<\gamma<1$ being CD; $0<\alpha<1,\gamma=1$ being CC. The stretched exponential relaxation pattern, known also as the Kohlrausch-Williams-Watts (KWW) function, has a simple form in time domain (1970) [29], namely it reads

[TABLE]

with $0<\alpha<1$. Here $\tau_{p}=\omega_{p}^{-1}$ is the time constant characteristic for a given material. The KWW function takes the simple exponential form if $\alpha=1$.

1.2 Experimental peculiarity of relaxation

In his two monographs A. K. Jonscher [10, 11] has shown that a common property of the empirical relaxation laws is that they exhibit the fractional-power dependence in the complex dielectric susceptibility $\chi(\omega)=\chi^{\prime}(\omega)-i\chi^{\prime\prime}(\omega)$, i. e.

[TABLE]

and

[TABLE]

where positive constant $\omega_{p}$ is the loss peak frequency characteristic for the investigated material, and $\Delta\chi(\omega)=\chi_{0}-\chi(\omega)$. It is worth noticing that this unique property is independent on any particular details of the examined systems and a great majority of dielectric materials may be characterized by two parameters (power-law exponents) $m$ and $1-n$, both falling strictly within the range (0,1]. This evidently shows Figure 1, where the exponents were defined for one hundred different materials (see Table 5.1 in [10]). In this connection it should be noticed that any satisfactory theory of relaxation must be capable of explaining this very general feature being so largely independent of the detailed physical nature of the materials involved.

Experimental data of relaxation can be characterized in both time and frequency domains. The inverse Stieltjes-Fourier transform

[TABLE]

relates the time-domain relaxation function $\phi(t)$ to the complex susceptibility $\chi(\omega)$ by the formula $\chi(\omega)=\phi^{*}(\omega)(\chi_{0}-\chi_{\infty})+\chi_{\infty}$, where $\phi^{*}(\omega)$ is the frequency domain relaxation function (shape function). The constant $\chi_{\infty}$ represents the asymptotic value of the susceptibility $\chi(\omega)$ at high frequencies, and $\chi_{0}$ is the value of the opposite limit. Note that the derivative $f(t)=-d\phi(t)/dt$, called the response function, is connected with the shape function $\phi^{*}(\omega)$ by the following relation

[TABLE]

that will be very useful for finding the latter function hereinafter. The existence of the following asymptotic responses, corresponding to Eqs.(3) and (4),

[TABLE]

in relaxation dynamics of complex systems has been established [10, 11, 30], the function $V(t)$ being regarded as

[TABLE]

The typical, fractional two-power-law behavior (3) and (4) is usually fitted with the HN function. In this case one has $m=\alpha$, $n=1-\alpha\gamma$ and $m\geq 1-n$. The exponents $m=1$ and $n=1-\gamma$ correspond to the CD relaxation, characterized by the short-time fractional power law only. The same property, however of different origins, concerns the KWW response for which the short-time power-law exponent $n=1-\alpha$.

Observe that fitting with the HN function of the atypical relaxation data (see Figure 1), for which the power-law exponents fulfill the opposite relation $m<1-n$, requires for the parameter $\gamma$ values greater than 1 [30]. As it will be shown below, the stochastic scenarios of relaxation do not allow us to derive this function with $\gamma>1$.

2 Definitions and terminology

2.1 Limit theorems

The probability theory considers a chance of the occurrence of an event in multiple repeating random experiments so as, for example, in a series of throws of a coin, where we can observe either its head or tail many times [31]. In stochastic modeling of kinetic processes the basic notation involves random variables. They are characterized by the distribution function $F(x)=\Pr(X\leq x)$ providing information about the probability $\Pr(x\leq X<x+dx)$, that the random variable $X$ takes a value between $x$ and $x+dx$, is equal to the difference $F(x+dx)-F(x)$. If the distribution function $F(x)$, $F(-\infty)=0$ and $F(+\infty)=1$, of the random variable $X$ fulfills the condition

[TABLE]

for every real $x$, then the function $f(x)$ is called the probability density function (pdf). The $n$th moment of a pdf $f(x)$ is the expected value of $X^{n}$, namely

[TABLE]

More generally, for any integrable function $g(\cdot)$ the expected value of $g(X)$ reads

[TABLE]

Note, the last definition will be often used in the present paper.

The study of sequences of random independent and identically distributed (iid) variables $X_{1},\,X_{2},\dots,\,X_{n}$ is one of the central topics in the probability theory. This is explained by some causes. At first, the statistical properties of the above sequence can be analyzed only asymptotically, i. e. when the number of variables $n\to\infty$. The distribution characteristics such as moments are calculated in this way. On the other hand, often the set $X_{1}$, $X_{2},\ \dots$ is a sequence of observations where the variable $X$ is observed repeatedly in time. Each individual observation is unpredictable, but the frequency of different outcomes over a large number of such observations becomes predictable. In particular, following the Bernoulli law of large numbers, in the experiments with only two results (“success” and “failure”) the frequency of the success will oscillate around the probability $p$ of the success [31]. The number of the success in $n$ trials is defined by the sum $\zeta_{n}=X_{1}+\dots+X_{n}$ having the binomial distribution. The strong law of large numbers states that the random variable $\frac{\zeta_{n}}{n}$ loses its randomness as the number $n$ of trials tends to infinity. Further studies of the deviation estimate $\left|\frac{\zeta_{n}}{n}-p\right|$ led to the first central limit theorem, i. e. the sum $\zeta_{n}$ for sufficiently large $n$, independently of distribution of a single component $X_{i}$ (but with the finite second moment), has the law close to the normal one. Namely, the central limit theorem answers why in so many uses (like the theory of errors, for example) one can find the probability distributions closely connected with the Gaussian one. Moreover, a wide circle of practical applications extends an essence of this theorem so that it has been generalized in many different ways. One of such generalizations concerns those distributions of $X_{i}$ that have no finite variance nor even mean value. Other direction to new limit theorems considers the operations on the sequence $X_{1}$, $X_{2},\dots$ different from the summation, as for example, in the extreme value theory [32] where the minimum and maximum operations are taken into account. Any case of the limit theorems indicates an asymptotic tendency of the sequence of random variables, as a result of an operation, to some non-degenerate random variable belonging to the class of limiting distributions (domain of attraction) different for every operation.

2.2 Lévy $\alpha$-stable distributions

In this subsection we present some basic facts on the Lévy $\alpha$-stable ($L\alpha S$ in short notation) distributions useful for the purpose of this article. The principal feature of these distributions is that they are completely described as limits of the normalized sums of iid summands [31]. Consequently, $L\alpha S$ distributions represent some kind of a universal law.

The distribution function $F(x)=\Pr(X\leq x)$ is called stable if for every $a_{1}>0$, $b_{1}$, $a_{2}>0$, $b_{2}$ there are constants $a>0$ and $b$ such that the equation

[TABLE]

holds. The symbol $F_{1}\ast F_{2}$ indicates the convolution of two distributions in the sense

[TABLE]

It turns out that always

[TABLE]

and the constant $\alpha$ is called the characteristic exponent of $L\alpha S$ distribution. Equation (8) can be solved in terms of characteristic functions, i. e., via Fourier transform

[TABLE]

For the distribution function $F(x)$ to be $L\alpha S$ it is necessary and sufficient that its characteristic function $h(s)$ is represented by the formula

[TABLE]

where $\alpha$, $\beta$, $\gamma$ and $\sigma$ are real constants with $\sigma\geq 0$, $0<\alpha\leq 2$ and $|\beta|\leq 1$. Here, $\alpha$ is the characteristic exponent, $\gamma$ and $\sigma$ determine location and scale. The coefficient $\beta$ indicates whether the $L\alpha S$ distribution is symmetric ($0<\alpha\leq 2$, $\beta=0$) or completely asymmetric ($0<\alpha<1$, $|\beta|=1$). The values $\alpha=2$ and $\beta=0$ yield the Gaussian distribution. As $h(s)$ is absolutely integrable, the corresponding $L\alpha S$ distribution has a density $f(x)$. Beautiful animations of the stable pdf’s with different values of the parameters are available on J. P. Nolan’s website (http://fs2.american.edu/jpnolan/www/stable/stable.html).

The most convenient formulation of the limit theorem, which gives description of the distribution law $F(x)$ governing the sum of a large number of mutually iid random quantities $X_{i}$, $i=1,2,\dots,n$, can be given in the following form : only L$\alpha$S distributions have a domain of attraction, i. e., there exist normalizing constants $a_{n}>0$, $b_{n}$ such that the distribution of $a_{n}^{-1}(X_{1}+X_{2}+\dots+X_{n})-b_{n}$ tends to $F(x)$ as $n\to\infty$. The normalizing constants can be chosen in such a way that $a_{n}=n^{1/\alpha}$.

Note, the random variable $X$ can describe an arbitrary physical magnitude (e. g. time, space, temperature, energy, etc). In particular, when $X$ is a waiting or residence time, the tail $\Pr(X>x)$ of the distribution $F(x)$ determines the survival probability. Let us add that the distribution of a non-negative random variable, say $X$, has a power asymptotic form if the tail $1-F(x)=\Pr(X>x)$ satisfies the condition

[TABLE]

for some $a>0$ and $\sigma>0$; that is, if for large values of $x$ the tail decays as a fractional power law $\sigma x^{-a}$. There are many different continuous and discrete distributions satisfying condition (11). Classical examples of continuous ones are the $L\alpha S$ laws, also the Pareto and Burr distributions with an appropriate choice of their parameters [31, 33].

If the distribution of random variable $X$ has a heavy tail with the parameter $0<a<1$, then the expected value $\langle X\rangle=\int x\,f(x)\,dx=\int x\,dF(x)$ is infinite. Note that in general if $\Pr(X>x)\sim 1/x^{a}$ for $x\to\infty$, then the moments $\langle X^{n}\rangle$ are finite for $a>n$. Therefore, the two considered attributes, the finiteness of the expected value and the heavy-tail property (11), clearly exclude each other. Besides, both provide only limited information on the corresponding distributions. Hence, the conditions put on the distributions of the microscopic quantities in the proposed scheme are rather general. On the other hand, by utilizing the limit theorems of probability theory the macroscopic result is determined in any detail.

2.3 Mixtures of distributions

Mixtures of distributions occur frequently in applications of the probability theory [31]. They also are directly relevant to problems of non-exponential relaxation. In this instance we deal with random variables the distribution of which depends on various factors, and all relaxing systems consist of many subsystems interacting among each other in a random way. Therefore, we call the sort of systems as a complex one. If $X$ is the random variable with pdf $f_{X}(x)$, then the random variable $aX$ ($a$ is constant) has the pdf $f_{X}(x/a)\frac{1}{a}$, and the random variable $X+a$ obeys the pdf $f_{X}(x-a)$. Let $Y$ be another random variable with pdf $f_{Y}(y)$. Now the product of random variables $XY$ takes the pdf in the integral form

[TABLE]

On the other hand, the pdf of the random variable $X/Y$ is written as

[TABLE]

The sum $X+Y$ is described by the convolution of pdfs, namely

[TABLE]

The relaxation rates of complex systems can depend on many parameters: temperature, defects, pressure and so on. Each of them has a very different distribution during a specific experimental scenario. However, macroscopic behavior of such systems is only a result of averaging such random effects. Thus, the mixtures of distributions become very helpful for the study of relaxation mechanisms.

2.4 Stochastic processes. Subordination

As J. L. Doob has defined [34], a stochastic process is “the mathematical abstraction of an empirical process whose development is governed by probabilistic laws”. There are two equivalent points of view about what is the stochastic process: (i) an infinite collection of random variables indexed by an integer or a real number often interpreted as time, and (ii) a random function of two or several deterministic arguments, one of which is the time $t$. It is convenient to consider separately the cases of discrete and continuous time. A discrete stochastic process $X=\{X_{n},\ n=0,1,2,\dots\}$ is a countable collection of random variables indexed by the non-negative integers, and a continuous stochastic process $X=\{X_{t},\ 0\leq t<\infty\}$ is an uncountable collection of random variables indexed by the non-negative real numbers. The Bernoulli process is perhaps the simplest non-trivial stochastic process. It is a sequence, $X_{0},X_{1},\dots$, of iid binary random variables that take only two values, 0 and 1. The common interpretations of the values $X_{i}$ are true or false, success or failure, arrival or no arrival, yes or no, etc. Note that the simple model of Bernoulli process initiated a great development of the studies on the limit theorems and served as the building block for other more complicated stochastic processes (Poisson process, renewal processes and others). The most known continuous stochastic process is the Brownian motion. Starting with 1827, when the botanist R. Brown observed zigzag, irregular patterns in the movement of microscopic pollen grains suspended in water, the phenomenon has found a satisfactory explanation only in 1905-1906 due to the physicists A. Einstein and M. Smoluchowski. Their probabilistic models have been based on assumption that the Brownian motion is a result of continual collisions of the pollen grains by the molecules of the surrounding water [35]. In 1923, the mathematician N. Wiener has proved the mathematical existence of Brownian motion as a stochastic process with the given properties. Any Brownian motion is a continuous time series of random variables whose increments are iid normally distributed with zero mean. It is plausible by the central limit theorem. Notice that this stochastic process is also a continuous-time analog to the simple symmetric random walk [35]. If one considers a massive Brownian particle under the influence of friction, the Ornstein-Uhlenbeck process has a bounded variance and admits a stationary probability distribution [36]. Eventually, this list of continuous stochastic processes unbarred doors to their study in different ways and under various conditions [37].

On the other hand, the diversity of stochastic processes may be extended notedly, if the parameter (index) $t$ varies stochastically. This approach, introduced by S. Bochner in 1949 [38] is called subordination [31, 39]. Then the process $Y(G(t))$ is obtained by randomizing the time variable of the stochastic process $Y(\tau)$ using a new “timer”, which is a stochastic process $G(t)$ with nonnegative independent increments. The resulting process $Y(G(t))$ is said to be subordinated to $Y(\tau)$ and is directed by the process $G(t)$, which is called a directional process. The directional process is often referred to as randomized time or operational time [31]. In general, the subordinated process $Y(G(t))$ can be non-Markovian, even if its parent process $Y(\tau)$ is Markovian.

3 Relaxation function as the initial-state survival probability.

Probabilistic models

The classical Debye relaxation law

[TABLE]

is characterized by a single relevant relaxation time $\tau_{p}$. To account for the non-exponential relaxation phenomenon, in the historically oldest attempt E. von Schweidler [18] assumed different parts of the orientational polarization to decline exponentially with different relaxation times $\tau_{i}$, yielding

[TABLE]

where the weights $p_{i}$ of the exponential decays fulfill

[TABLE]

Few years later, K.W. Wagner [40] proposed to use of a continuous distribution $w(\tau)$ of relaxation times

[TABLE]

where $\int_{0}^{\infty}w(\tau)\,d\tau=1$.

This approach is microscopically arbitrary since it does not yield any constraints on the random microscopic scenario of relaxation. The probability density function $w(\tau)$, as well as, the weights $p_{i}$ are determined only by the empirical patterns of $\phi(t)$. This simple way to derive the non-exponential decay is associated with a picture of parallel relaxations, in which each degree of freedom (each relaxation channel) relaxes independently with random relaxation time [42, 43, 41, 44, 45, 46, 47, 6, 48, 49, 7, 50]. From the probabilistic point of view, both above formulas reveal the weighted average $\left\langle\dots\right\rangle$ of an exponential relaxation

[TABLE]

with respect to the distribution of the random effective relaxation time $\widetilde{T}$ with support of $\tau\in\left[0,\infty\right)$.

Contrary to models that were based on a parallel addition of relaxation contributions, the model presented in [51] proposes a serial summation of a hierarchy of relaxations extending over the same spatial range. The authors pointed out that a group of dipoles must adopt to a specific configuration before a subset can relax, which then releases the constraints preventing a further subset from relaxing, and so on. Although it has been realized in many approaches that the individual dipoles and their environment do not remain independent during the regression of fluctuation, as yet no microscopic model has been based directly on this conclusion. The exception is the cluster model [52, 2, 53, 54, 55, 56], which derived entirely new expressions from a consideration of the way in which the energy contained in fluctuation is distributed over a system of interacting clusters. This is also the only theory in which the results obtained are in agreement with empirical functions input to fit the experimental data for $\phi(t)$ in the short- $t\ll\tau_{p}$ and the long-time $t\gg\tau_{p}$ limits.

3.1 Transition probability

Relaxation properties of dielectric materials have been the subject of experimental and theoretical investigations for many years. This is not only due to the need for understanding of the electrical properties of various technological materials, but it has also been realized that the basic physics of the dielectric relaxation response leads to interesting questions about the theoretical description of physical relaxation mechanisms in disordered (complex) systems.

Relaxation phenomena are experimentally observed when a physical macroscopic magnitude (concentration, current, etc.), characteristic for the investigated system, monotonically decays or grows in time. In case of the dielectric relaxation this process is commonly defined as an approach to equilibrium of a dipolar system driven out of equilibrium by a step or alternating external electric field. The time-dependent response of relaxing systems to a steady electric field is described by the relaxation function $\phi(t)$ (satisfying $\phi(0)=1$ and $\phi(\infty)=0$) which is a solution of the two-state master equation

[TABLE]

The non-negative quantity $r(t)$ is the transition rate of the system (i. e. the probability of transition per unit time), see e. g. [42]. Consequently, the function $\phi(t)$ has a meaning of the survival probability $\Pr(\tilde{\theta}\geq t)$ of a non-equilibrium initial state of the relaxing system [7]. In other words, $\phi(t)$ is determined by the probability that the system as a whole will not make a transition out of its original state for at least time $t$ after entering it at $t=0$.

Note, if instead of the decay of polarization one observes its increase under the influence of the steady external electric field, the proper function to describe this relaxation process would be just $F(t)=1-\phi(t)$ satisfying $F(0)=0$ and $F(\infty)=1$ (see, e.g. [10]). From the probabilistic point of view, the response function $f(t)$ is just the probability density function of the waiting-time distribution $F(t)=1-\phi(t)=\Pr(\tilde{\theta}<t)$.

Let us assume that a relaxing physical system undergoes an irreversible transition from initial state ${\bf 1}$, imposed at time $t=0$, to state ${\bf 2}$ that differs from ${\bf 1}$ in some physical parameter. The transition ${\bf 1}\to{\bf 2}$, defined as a change of this parameter, takes place at random instant of time. To be precise, we consider the conditional probability $p(t,dt)$ that the system as a whole will undergo the transition during a time interval $(t,t+dt)$ provided that the transition did not occur before time $t$, i. e.

[TABLE]

where $\tilde{\theta}$ is the system’s random waiting time for transition ${\bf 1}\to{\bf 2}$. The conditional probability $p(t,dt)$ defined in (16) can be expressed as

[TABLE]

where $\Pr(\tilde{\theta}\geq t+dt)$ and $\Pr(\tilde{\theta}\geq t)$ are the survival probabilities, i. e. the probabilities that the system will remain in state ${\bf 1}$ until time $t+dt$ and $t$, respectively. As $dt\to 0$, probability (17) can be rewritten in a form more useful for further considerations, namely

[TABLE]

In this case the function

[TABLE]

is nothing else than the intensity of time-dependent transition probability (transition rate) [57, 58]. Then the survival probability $\Pr(\tilde{\theta}\geq t)=\phi(t)$ can be found at once if one knows the explicit form of the intensity, namely

[TABLE]

As the intensity of transition is essentially time dependent, the evolutionary law for the entity is of a non-exponential form. For the time-independent intensity $r(s)=\tau_{p}^{-1}={\rm const}$ one obtains

[TABLE]

what recovers the classical exponential decay law (12) with the value $\tau_{p}^{-1}$ determining the relaxation rate of the transition process.

Unfortunately, the time-dependent relaxation rate $r(t)$ in Eq. (15) can take the very cumbersome form for empirical laws of relaxation (especially, for the HN one). Moreover, the determinism appears from empirical laws, observed on the macroscopic level, while randomness is induced by variations in the local environment. In order to understand the phenomenon of the universal relaxation response one needs to consider relaxation in complex systems in a way that separates it from a particular physical context, and observation that the dynamics of such systems is characterized by seemingly contrary states, i. e., local randomness and global determinism, is crucially relevant to this issue. These two states, in a natural way, can coexist in the framework of the limit theorems of probability theory when the relaxation function $\phi(t)$, being the survival probability of the system in the initially imposed non-equilibrium state until time $t$, is described by the first-passage of the system.

3.2 Microscopic stochastic scenario

In any complex system, capable of responding to an external electric field, it is possible that only a part of the total number $N$ of dipoles in the system is able to follow significantly changes of the field [10]. In general, the distribution of the random waiting time $\tilde{\theta}$ of the entire system is determined by the first passage of the system from its initial state [59, 58].

Consider (as in the preceding subsection) a system of $N$ entities, each waiting for transition ${\bf 1}\to{\bf 2}$ for a random time $\theta_{i}$, where $1\leq i\leq N$. Generally speaking, the waiting times $\theta_{1},\theta_{2},\ldots,\theta_{N}$ form an arbitrary sequence of iid non-negative random variables. The entities undergo transition in a certain order that can be reflected in the notion of order statistics, that provide a non-decreasing rearrangement $\theta_{(1)}\leq\theta_{(2)}\leq\ldots\leq\theta_{(N)}$ of times $\theta_{1},\theta_{2},\ldots,\theta_{N}$ [58, 31]. From this rearrangement follow two obvious statistics, the first and the $N$th one, $\theta_{(1)}$ and $\theta_{(N)}$, respectively. Now, denoting an unknown (random) number of individual transitions occurred until time $t>0$ by $\eta_{N}(t)$, we can connect the event number $\{\eta_{N}(t)=n\}$ with the order statistics via the relations

[TABLE]

The first of them indicates that no transition has occurred in this system until time $t$. The second one shows a tendency to ${\bf 1}\to{\bf 2}$. The population of entities in the state ${\bf 1}$ are decreased step-by-step in favor of the relaxation output ${\bf 2}$. The last expression of Eq.(19) means that all transitions have been finished up to time $t$.

Let us now introduce a notion of the initial-state survival probability of the entire $N$-dimensional system as {\rm Pr}\bigl{(}\widetilde{\theta}_{N}\geq t\bigr{)}, i. e. the probability that transition of the system as a whole has not occurred prior to a time instant $t$, where $\widetilde{\theta}_{N}$ denotes the system’s waiting time for transition from its initial, imposed state. The probability {\rm Pr}\bigl{(}\widetilde{\theta}_{N}\geq t\bigr{)} means also that there is no any individual transition until time $t$. Therefore, from Eq.(19) we have

[TABLE]

As a rule, the macroscopic systems consist of a large number of relaxing entities so that the relaxation function can be approximated by the weak limit in distribution

[TABLE]

where $A_{N}$ denotes a sequence of normalizing constants, and $\stackrel{{\scriptstyle d}}{{=}}$ means “equal in law”. In relation to the above definition of the relaxation function, the frequency-domain shape function (6) can be written as

[TABLE]

where $\langle\dots\rangle$ denotes an average with respect to the distribution of the system’s effective waiting time $\widetilde{\theta}$. It follows from limit theorems of the extremal value theory [32], that since the sequence of waiting times $\theta_{1},\theta_{2},\ldots,\theta_{N}$ consists of iid non-negative random variables, the above definition of the relaxation function leads to the result

[TABLE]

where $A$ and $\alpha$ are positive constants. Observe that this form of the relaxation function, being just the tail of the well-known Weibull distribution [31], contains three possible cases: the stretched exponential behavior if $0<\alpha<1$, the exponential if $\alpha=1$, and the compressed exponential one if $\alpha>1$. At this point it is natural to ask how, within the proposed stochastic scenario, one can derive the stretched exponential KWW function, as well as, the other empirical relaxation patterns (see Eq.(1)). To solve this problem one has to realize that, in general, in empirical relaxation evidence we observe two classes of relaxation responses. Namely, a class exhibiting the short-time power law only (fitted with the KWW or CD function), and a class exhibiting both short- and long-time power laws (fitted with the HN function). Hence, the first step to solve this problem is to find a rigorous mathematical condition yielding the KWW short-time power law in the framework of the above stochastic scenario.

3.3 Stretched exponential relaxation

Traditional interpretation of the non-exponential relaxation phenomena is based on the concept of a system of independent, exponentially relaxing species (dipoles) with different relaxation rates [41]. The exponential relaxation of an individual dipole in this case is conditioned only by the value taken by its relaxation rate. So, taking into account influence of the local random environment on the entity, one can conclude [59, 58]: if the relaxation rate $\beta_{i}$ of ith dipole has taken the value $b$, then the probability that this dipole has not changed its initial aligned position up to the moment $t$, is

[TABLE]

The random variable $\beta_{i}$ denotes the relaxation rate of ith dipole and the variable $\theta_{i}$, the time needed for changing its initial orientation; $\beta_{1},\beta_{2},\dots$ and $\theta_{1},\theta_{2},\dots$ form sequences of non-negative, iid random variables. The randomness of the individual relaxation rate is motivated by the fact that in a complex system its entity can be into many states or even pass through a whole hierarchy of substates within states, and the distribution of individual relaxation rates effectively accounts for the transition intensity between the states and substates.

From the law of total probability [31], we have

[TABLE]

where $F_{\beta_{i}}(b)$ is the distribution function of each relaxation rate $\beta_{i}$. In other words, $F_{\beta_{i}}(b)$ denotes the probability that the relaxation rate of ith dipole has taken a value less than or equal to $b$. Formula (24) shows that if one takes into account influence of the local random environment on relaxation behavior of a dipole, its initial-state survival probability decays non-exponentially. Only if the influence is deterministic, i.e., the individual relaxation rate $\beta_{i}$ takes the value $b_{0}$ with probability 1, given by the pdf of the Dirac $\delta$-function form $dF_{\beta_{i}}(b)=\delta(b-b_{0})\,db$, the individual survival probability decays exponentially $\Pr(\theta_{i}\geq t)=e^{-b_{0}t}$.

In order to obtain an explicit form of the relaxation function $\phi(t)$ defined in (21), let us observe that the right-hand expression in (24) is just the Laplace transform of the distribution function $F_{\beta_{i}}(b)$ at point $t$,

[TABLE]

Because $\theta_{i}$ are independent random variables, we get

[TABLE]

The Nth power of the Laplace transform of the non-degenerate distribution function converges to the non-degenerate limiting transform, as $N$ tends to infinity, if and only if $F_{\beta_{i}}(b)$ belongs to the domain of attraction of the completely asymmetric $L\alpha S$ law $F_{\widetilde{\beta}}(b)$. Then, for some $\alpha$, $0<\alpha<1$, we have

[TABLE]

where $A$ is a positive constant. Hence, the limiting transform in (25) is the Laplace transform of the $L\alpha S$ distribution with the non-negative support and the stable parameter $\alpha$ belonging to the range (0,1).

It is not necessary to know the detailed nature of $F_{\beta_{i}}(b)$ to obtain the above stretched exponential (KWW) limiting form. In fact, this is determined only by tail behavior of $F_{\beta_{i}}(b)$ for large $b$, see Eq.(11), and so a good deal may be said about the asymptotic properties based on rather limited knowledge of the properties of $F_{\beta_{i}}(b)$. In other words, the necessary and sufficient condition for the relaxation rate $\beta_{i}$ to have the limiting transform in (25) is the self-similar property in taking the value greater than $b$ and the value greater than $xb$, where $x$ is a positive constant, and $b$ takes a large value. It has been suggested [55, 7] that self-similarity (fractal behavior) is a fundamental feature of relaxation in real materials. This result, obtained here by means of pure probabilistic techniques, independently of the physical details of dipolar systems, is in agreement with models [7, 50, 51] identifying this region of fractal behavior.

Let us observe that the right-hand side of formula (24) can be also interpreted as the weighted average of an exponential decay

[TABLE]

where the mean value $\langle\dots\rangle$ is taken with respect to the relaxation-rate probability distribution $F_{\beta_{i}}(b)$. This leads to

[TABLE]

where $(\beta_{1},\beta_{2},...\beta_{N})$ are the non-negative iid random relaxation rates of individual transitions. If $\beta_{i}$ has a finite mean, i. e. $\langle\beta_{i}\rangle<\infty$, then the macroscopic development gives nothing new because the relaxation evolves exponentially with a constant rate $b_{0}$, and $\langle\widetilde{\beta}_{N}=\sum^{N}_{i=1}\beta_{i}/N\rangle=b_{0}$ as $N\to\infty$. But the stochastic picture changes drastically, if the sum

[TABLE]

consists of rates no having any mean. Summation of iid random variables is well known in literature [31, 60] and the resulting completely asymmetric $L\alpha S$ distribution $F_{\widetilde{\beta}}(b)$ of the effective relaxation rate $\tilde{\beta}$ can be approximated by the weak limit

[TABLE]

In practice, even $N\approx 10^{5}-10^{6}$ can suffice to replace adequately $\widetilde{\beta}_{N}$ in (26) by the limit $\widetilde{\beta}$. Taking into account Eqs. (20)-(28), we get

[TABLE]

what again yields the KWW stretched exponential decay (25).

Therefore, the relaxation function (21) with $A_{N}=N^{1/\alpha}$, for some $0<\alpha<1$, is well defined and equals

[TABLE]

where $A=\tau_{p}^{-1}$ (see Eq.(2)). When $\alpha\to 1$, the theoretically derived KWW function (30) obtains the D form (12). From the mathematical point of view [31, 60, 61] this corresponds to the case of degenerate distribution function $F_{\widetilde{\beta}}(b)$, i. e. to the case when the effective random relaxation rate $\widetilde{\beta}$ can take only one value. The corresponding pdf is then of the Dirac $\delta$-function form. At this point we have to stress that the degenerate distributions (of different, studied below physical magnitudes) yield the limiting value 1 of the HN and KWW exponents (see Eqs. (1) and (2)). So, to avoid confusion between the theoretical (0,1) and the experimental (0,1] ranges of possible values, taken by the characteristic exponents, we will always in our theoretical studies include the degenerate distributions.

Following the historically oldest approach to non-exponential relaxation [41], the relaxation function can be expressed as $\phi(t)=\langle e^{-t/\widetilde{T}}\rangle$, since it has been assumed that non-exponential relaxation function takes the form of a weighted average of an exponential decay $e^{-t/\tau}$ with respect to the distribution $w(\tau)\,d\tau$ of the random effective relaxation time $\widetilde{T}$, see Eq.(13). As the effective relaxation rate $\widetilde{\beta}=1/\widetilde{T}$, the formula (13) can be rewritten as follows

[TABLE]

where $\left(\widetilde{\beta}:b\in[0,\infty)\right)$. This representation assigns any non-exponential relaxation function $\phi(t)$ to the Laplace transform of the effective relaxation-rate distribution $g(b)$. The probability density functions $w(\tau)$ (see Eq.(13)) and $g(b)$ (see Eq.(31)) are related to each other, namely $g(b)=b^{-2}w(b^{-1})$. The relationship between $g(b)$ and $w(\tau)$, corresponding to the KWW relaxation, allows us to show [50] that in contrast to the momentless distribution $dF_{\widetilde{\beta}}(b)=g(b)\,db$ of the effective relaxation rate $\widetilde{\beta}$, the distribution $dF_{\widetilde{T}}(\tau)=w(\tau)\,d\tau$ possesses finite average and higher moments of effective relaxation time $\widetilde{T}$. Notice, the relaxation rates are additive, but the relaxation times are not. Therefore, the relaxation rates as random variables are more convenient for the probabilistic formalism based on the limit theorems of probability theory. Hence, in further study only formula (31) will be utilized.

Let us observe that independently on a statistical distribution of relaxation rates $\beta_{i}$ we find in expression (23) a hidden assumption. Namely, each relaxing dipole after a sufficiently long time (after removing the electric field) changes its initial position with probability 1, i. e.

[TABLE]

Such an assumption is the main reason why the relaxation function (21) cannot have any other form than the KWW one (30). The above analysis gives also an insight into the physical origins of the short-time power law observed in all non-exponential relaxation responses. For the simplest non-exponential case (30), the response function reads

[TABLE]

where $n=1-\alpha$ results from the $L\alpha S$ distribution of the effective relaxation rate $\widetilde{\beta}$. The power-law exponent $n$ is determined by the long-tailed properties (see Eq.(11)) of this distribution $1-F_{\widetilde{\beta}}(b)\propto b^{-\alpha}$ for $b\to\infty$.

In order to obtain a class of dielectric responses exhibiting both the short- and the long-time power law, one should modify either the assumption (23) to define the random waiting time which can be infinite with some non-zero probability or modify the definition of the relaxation function (21) to account for the random number of individual relaxation contributions (relaxation channels). As we will see below, such a modification, being in agreement with physical intuition on relaxation mechanisms, leads us directly to the non-exponential responses (7). In proposed schemes of relaxation the KWW and D functions are included as special cases.

3.4 Conditionally exponential decay model

Let us assume independent exponential relaxations constrained by the maximal time of a structural reorganization in all surrounding clusters (each consisting of a dipole and its non-polar environment). In a system composed of $N$ of relaxing dipoles, the probability [62, 63] that the ith dipole has not changed its initial position up to the moment $t$ equals $\exp[-b\min(t,s)]$, if its relaxation rate has taken the value $b$ and the maximal time of the structural reorganization $\eta_{i,\rm max}=a_{N}^{-1}\max(\eta_{1},\dots,\eta_{i-1},\eta_{i+1},\dots,\eta_{N})$ in all surrounding clusters (under the suitable normalization) has been equal to $s$, i. e.

[TABLE]

for $b>0,s>0,t\geq 0$. The random variable $\beta_{i}$ denotes the relaxation rate of the ith dipole and the variable $\eta_{i}$, the time needed for the structural reorganization of ith cluster. The variable $\theta_{iN}$ denotes the time needed for changing the orientation by the ith dipole in the system consisting of $N$ relaxing dipoles. $\beta_{1},\beta_{2},\dots$ and $\eta_{1},\eta_{2},\dots$ form independent sequences of non-negative, iid random variables. The variables $\theta_{1},\dots,\theta_{N}$ are also non-negative, iid for each $N$. It follows from (33) that the random variable $\theta_{i}$ depends on the random variable $\beta_{i}$ and on the sequence of random variables $\eta_{1},\dots,\eta_{i-1},\eta_{i+1},\dots,\eta_{N}$.

In contrast to (32) we have

[TABLE]

It means that dipoles altered by the external electric field do not have to change their initial positions with probability 1 after removing the field as $t$ tends to infinity (with some probability their initial states are “frozen”). In this case, because of the improper form of the distribution (33), the relaxation function (21) cannot be expressed as the waiting average with respect to the effective relaxation rate distribution (see Eq.(31)). Instead, a general relaxation equation, fulfilled by function (21), can be derived [63].

Since sequences $\beta_{1},\beta_{2},\dots$ and $\eta_{1},\eta_{2},\dots$ are independent, we have from the law of total probability

[TABLE]

where $F_{\eta,N}(s)$ denotes the distribution function of the random variable which has the form $a_{N}^{-1}\max(\eta_{1},\dots,\eta_{i-1},\eta_{i+1},\dots,\eta_{N})$, i. e. the probability that this random variable has taken a value less than or equal to $s$. Since $\eta_{j}$ are iid random variables, we have $F_{\eta,N}(s)=[F_{\eta}(a_{N}s)]^{N-1}$, where $F_{\eta}(s)$ denotes the distribution function of each $\eta_{j}$. Assuming $F_{\eta}(s)$ differentiable, we have $F_{\eta,N}(s)$ differentiable, too, and

[TABLE]

From the law of total probability once again, and from the Lebesgue theorem [31], we have

[TABLE]

where ${\cal L}(F_{\beta_{i}}(b);t)$ is the Laplace transform of the distribution function $F_{\beta_{i}}(b)$ of each $\beta_{i}$ at the point $t$.

Because $\theta_{iN}$ are iid random variables for each $N$, we have

[TABLE]

Using the mathematical trick

[TABLE]

it follows from (34) that

[TABLE]

As we know from the preceding section, the $N$th power of the Laplace transform of a non-degenerate distribution function $F_{\beta_{i}}(b)$ converges to the non-degenerate limiting transform, as $N$ tends to infinity, if and only if $F_{\beta_{i}}(b)$ belongs to the domain of attraction of the $L\alpha S$ law, and, for some $0<\alpha<1$, we have

[TABLE]

where $A$ is a positive constant. At the same time, the value

[TABLE]

tends to a non-degenerate distribution function of non-negative random variable, as $N$ tends to infinity, if and only if $F_{\eta_{i}}(s)$, the distribution function of each $\eta_{i}$, belongs to the domain of attraction of the max-stable law of type II [32]. Then, for the normalizing constant $a_{N}$ proportional to $N^{1/\alpha}\inf\{t:F_{\eta}(t)\geq 1-(1/N-1)\}$ we have

[TABLE]

for some positive constants $\alpha$ and $\kappa$, and $A$ taken from Eq.(37). To obtain the limiting forms (37) and (38) we need not know the detailed nature of $F_{\beta_{i}}(b)$ and $F_{\eta_{i}}(s)$. In fact, this is determined only by the behavior of the tail of $F_{\beta_{i}}(b)$ for large $b$ and of the tail of $F_{\eta_{i}}(s)$ for large $s$, i. e. the necessary and sufficient conditions for the relaxation rate $\beta_{i}$ and for the structural reorganization time $\eta_{i}$ to have the limits in Eqs. (37) and (38) are the self-similar properties, firstly of $\beta_{i}$, in taking the value greater than $b$ and the value greater than $xb$, and secondly of $\eta_{i}$ in taking the value greater than $s$ and the value greater than $xs$.

The relaxation function in Eq.(21) with $A_{N}=N^{1/\alpha}$ is well defined and, by Eqs. (35)–(38), fulfills the general relaxation equation (a kinetic equation with a time-dependent transition rate $r(t)$, see Eq.(15))

[TABLE]

Recall that the parameter $A$ has the sense of $\tau_{p}^{-1}$. The coefficient $\kappa$ is a consequence of normalization in the limiting procedure in Eq.(38). It decides how fast the structural reorganization of clusters is spread out in a system; $\kappa\to 0$ means the case in which cluster structure is neglected. If $\kappa\to 0$, Eq.(39) takes the well-known form [55, 7, 42]

[TABLE]

with the solution (30). In the general case we get the solution in an integral form

[TABLE]

where $c=\kappa^{-1}$ and

[TABLE]

A similar form has been obtained as a result of the studies of different approaches (the Förster direct-transfer model, the hierarchically constrained dynamics model, and the defect-diffusion model) analyzing non-exponential relaxations, with emphasis on the stretched exponential KWW form [7, 64, 65]. Although each model describes a different mechanism, they have the same underlying reason for the stretched exponential pattern: the existence of scale invariant relaxation rates. Presenting one more approach, we have obtained the KWW relaxation function (30) as a special case of Eq.(39) when $\kappa\to 0$. We have also shown that the underlying reason for this is the existence of a type of self-similarity in the behavior of relaxation rates.

For practical purposes, according to [66], the solution of Eq.(39) can be presented in the following form

[TABLE]

where $\Gamma(a,z)$ is the incomplete gamma function defined as

[TABLE]

It follows from Eq.(39) that the relaxation response may be written as

[TABLE]

Then, for the short-time regime its asymptotic behavior is

[TABLE]

since $\phi(0)=1$. On the other hand, the long-time trend follows

[TABLE]

where $\gamma_{E}\approx 0.577216$ is the Euler constant [66]. Thus, the response function $f(t)$ can exhibit the power-law properties in both short- and long-time limits, namely

[TABLE]

where $n=1-\alpha$ and $m=\alpha/\kappa$. The relaxation function $\phi(t)$ is determined by three parameters: $0<\alpha<1$, $A>0$ and $\kappa>0$. The parameter $\kappa$ distinguishes the fractional two-power-law behavior from the one-power-law KWW response, i. e. if $\kappa$ is small, the general relaxation solution of Eq.(39) takes the form which is just the KWW relaxation function. Moreover, if $\alpha\to 1$ we obtain the rarely observed D case. For $\kappa<1$ the relaxation function of Eq.(39) describes the typical case of Figure 1, and for $\kappa>1$ we obtain the less typical relaxation behavior. Note also that the formalism of coupled cluster interactions finds a good support in experimental studies [67, 68, 69].

3.5 Relaxation of hierarchically clustered systems

The cluster model concept [2, 53] presents a radical departure from the traditional interpretation of relaxation based on independent exponentially relaxing entities. The realistic idea originates from imperfectly ordered states of complex systems and their evolution. In this case the systems, which exhibit position or orientation relaxation, are composed of spatially limited regions (clusters). Because the structural order within any cluster is incomplete, there are internal and external dynamics of clusters. When an external field acts on such a system, entities of this system take positions along the field direction, but the positions will be very dependent on the local structure of the system, i. e. from defects of different types. With regard to the imperfect structure the arrangement of entities after removing the external field starts to lose spatial uniformity. During this process of relaxation to an equilibrium geometry, the strongly coupled local motions are expected to arise firstly, thereby breaking down the arrangements into clusters that leads to weakly coupled inter-cluster motions forming a constraint hierarchy of interacting clusters and their long-range compositions. Each of these processes have their own characteristic contribution to the macroscopic evolution of the system as a whole. This is the reasonably natural way to modify a traditional approach of independent relaxing entities with random relaxation rates into a multilevel summation of a hierarchy of cluster relaxations with their random relaxation rates.

3.5.1 Havriliak-Negami function

Before going into details of the random-cluster relaxation model [70, 71, 72] let us first discuss the stochastic representation (22) of the HN function. Take into account the random effective waiting time $\widetilde{\theta}$ for transition of the relaxing system [73] as a mixture of random variables

[TABLE]

Here $S_{\alpha}$ is such a positive random variable that its Laplace transform is the stretched exponential function

[TABLE]

with $0<\alpha<1$. It is a well-known fact [31] that in the above relation the random variable $S_{\alpha}$ has to be distributed according to the completely asymmetric $L\alpha S$ law with the pdf $h_{\alpha}(t)$ (for details see [60, 61]). The pdf of the random variable $S_{\alpha}$ tends to the degenerate form $h_{1}(t)=\delta(t-1)$ (given by the Dirac $delta$-function) as $\alpha\to 1$.

The positive random variable $\varGamma_{\gamma}$ in Eq.(45) is independent of $S_{\alpha}$ and distributed according to the gamma law $G_{\gamma}(t)$ defined by the pdf of the form

[TABLE]

with $\Gamma(\gamma)$ being the Euler’s gamma function [33]. It is worth noting that the Laplace transform of $\varGamma_{\gamma}$ reads

[TABLE]

Using the properties (46) and (47), we derive the explicit form of Eq.(22)

[TABLE]

Therefore, the waiting time $\widetilde{\theta}_{\rm HN}$ with $0<\alpha,\gamma<1$ represents the HN relaxation pattern and, moreover, in the time domain we have

[TABLE]

where $G_{\gamma}(x)=\int_{0}^{x}g_{\gamma}(t)\,dt=1-\Gamma(\gamma,x)/\Gamma(\gamma)$, and $\Gamma(\gamma,x)$ is the upper incomplete gamma function.

On the other hand, the distribution of the random variable $\widetilde{\theta}_{\rm HN}$ can be identified as the Mittag-Leffler distribution [74] and, hence, the time-domain HN relaxation function is represented by the following series

[TABLE]

where $E^{\gamma}_{\alpha,\beta}(x)$ is the three-parameter Mittag-Leffler function [75].

Formula (45) and hence the time-domain relaxation function (48) take on simpler forms in case of the CD, CC and D responses. For the CD function ($\alpha=1$, $0<\gamma<1$) one gets the gamma waiting-time distribution and

[TABLE]

The respective response function equals $f_{\rm CD}(t)=g_{\gamma}(t/\tau_{p})/\tau_{p}$. Similarly, in case of the D function ($\alpha=\gamma=1$) the corresponding waiting time is distributed according to the exponential law and

[TABLE]

For the CC response ($0<\alpha<1$, $\gamma=1$) the gamma random variable becomes an exponential one $\varGamma_{\gamma}=\varGamma_{1}$ and the series representation (49) simplifies to the one-parameter Mittag-Leffler function

[TABLE]

It is worth noting that all the waiting-time distributions, underlying the considered empirical relaxation responses, are infinitely divisible [31, 73].

3.5.2 Infinite mean cluster sizes

Let us now study the complex dynamics of clustered systems from the probabilistic point of view [70, 71, 72]. In every complex system capable of responding to an external field, the total number $N$ of entities in the system may be divided into two parts. One of them includes so-called active entities being able to follow changes of the field. Another part consists of inactive neighbors. Even if some entities do not contribute directly to the relaxation dynamics, they may affect the stochastic transition of the active ones. Note, this influence lies in the properties of individual relaxation rates $\beta_{1N},\beta_{2N},\dots$ of the active entities in the system. According to the concept of relaxation rates, the individual rates take the form $\beta_{iN}=\beta_{i}/A_{N}$, where $\beta_{i}$ is independent of $N$, and $A_{N}$ is the same normalizing constant for each entity. Assume further that the $i$th active entity interacts with $N_{i}-1$ inactive neighbors forming a cluster of size $N_{i}$. The unknown number $K_{N}$ of active entities in the system, random in general, equals to the number of clusters due to the local interactions. The value $K_{N}$ is determined by the first index $k$ for which the sum $N_{1}+\dots+N_{k}$ of the cluster sizes exceeds $N$, the total number of entities. Mathematically this proposition can be written as

[TABLE]

where $\{k:X\}$ implies the value of $k$ such that $X$ holds. Interactions among active entities have a local character because of their surrounding by inactive entities (due to screening effects, for example, [10]). Therefore, every evolving active entity may “feel” only some of other active neighbors. In these conditions nothing prevents the emergence of cooperative regions (super-clusters) built up from the active entities and their surroundings. Let the random number $L_{N}$ of such super-clusters be determined by their sizes $M_{1},M_{2},\dots$. In the same way as (51) we define

[TABLE]

where $M_{j}$ is a number of interacting active entities in the j th super-cluster. A contribution of each super-cluster to the total relaxation rate is the sum of the contributions of all active entities in the super-cluster. Hence, for the j th super-cluster, its relaxation rate, say $\overline{\beta}_{jN}$, is equal to

[TABLE]

For $j=1$ it is simply the sum

[TABLE]

for $j=2$ it becomes

[TABLE]

and so on. The effective representation of the system as a whole is provided by the total relaxation rate $\widetilde{\beta}_{N}$ which is the sum of the contributions over all super-clusters

[TABLE]

In fact, considering relaxation phenomena, one usually deals with systems consisting of a large number of relaxing entities so that the weak limit

[TABLE]

can describe the entire system, and, hence the relaxation function reads

[TABLE]

In general, all the quantities $N_{i}$, $M_{j}$, $\beta_{iN}$, together with those defined by them, must be considered as random variables. The point is that the number of relaxing entities directly engaged in the relaxation process, their locations, as well as their “birth” and “death”, are random. Obviously, their stochastic features would determine the properties of the total relaxation rate $\widetilde{\beta}_{N}$ if they were known. But they are rather not known. Nevertheless, on the basis of the limit theorems of probability theory, it is possible to define the distribution of the limit $\widetilde{\beta}$ representing a macroscopic relaxing system, even with rather limited knowledge about the distributions of micro/mesoscopic random quantities used in the model.

Assume that $M=\{M_{j},\ j=1,2,\dots\}$, $N=\{N_{i},\ i=1,2,\dots\}$, and $\beta=\{\beta_{i},\ i=1,2,\dots\}$ are independent sequences, each of which consists of iid positive random variables, $M_{j}$ and $N_{i}$ being integer-valued. For the sake of simplicity, let us introduce the following notations

[TABLE]

for $k=1,2,\dots$, and $X$ corresponding to one of sequences $M$, $N$, or $\beta$. Denoting the sum of the contributions over all super-clusters (56) for convenience as $\widetilde{\beta}_{n}$, the latter takes the form of a normalized random sum

[TABLE]

with the random index $\Sigma_{n}={\cal S}_{M}(\nu_{M}(\nu_{N}(n)))$ independent of the components $\beta_{j}$, and $a_{n}$ as a sequence of the normalizing constants. If the distribution of a positive random variable $X_{j}$ has a heavy tail as defined in (11), then asymptotic properties of $\widetilde{\beta}_{n}$ (as $n\to\infty$) can be found exactly.

Assume that both $N_{j}$ and $\beta_{j}$ have heavy-tailed distributions with the same exponent $\alpha\in(0,1]$. Following [31], we have

[TABLE]

where random variables $S_{\alpha}$ and $Z_{\alpha}$ are distributed with completely asymmetric $L\alpha S$ distributions, $c_{1}$ and $c_{2}$ being positive constants. Here $\stackrel{{\scriptstyle d}}{{\rightarrow}}$ means “tends in distribution”. Since $\nu_{N}(n)$ and ${\cal S}_{N}(k)$ are connected by the relation

[TABLE]

for any $n,k=1,2,\dots$, then

[TABLE]

If the distribution of $M_{j}$ has also a heavy tail, according to Eq.(11), with $\gamma\in(0,1]$, as has been shown in [31], the normalized random sum

[TABLE]

for $n\to\infty$ behaves as a reciprocal to the beta-distributed random variable ${\cal B}_{\gamma}$ governed by the following pdf

[TABLE]

where $\Gamma(\cdot)$ is the gamma function. It should be pointed out that the distribution is a particular case of the beta-distribution [33]. As the random sequences $\{{\cal S}_{M}(\nu_{M}(n))\}$ and $\{\nu_{N}(n))\}$ are independent, it follows from [76] that the results (59) and (60) allow us to get

[TABLE]

Using the relations (58) and (61), we find the tendency of the normalized random sum $\widetilde{\beta}_{n}$ for $n\to\infty$, namely

[TABLE]

Derivation of $\lim_{n\to\infty}\widetilde{\beta}_{n}$ has been presented in [77].

Recall, the relaxation response can be associated with $\widetilde{\theta}$, the system’s waiting time for its transition from the initially imposed state, and $\widetilde{\beta}$, the effective relaxation rate. The random variables $\widetilde{\theta}$ and $\widetilde{\beta}$ are strictly connected with each other, as $\phi(t)=\Pr(\widetilde{\theta}\geq t)={\cal L}(\widetilde{\beta};t)$, where ${\cal L}(X;t)=\left\langle\exp(-Xt)\right\rangle$ denotes the Laplace transform of the pdf of a random variable $X$. As it has been shown above (see Eqs.(45)-(49)), the relaxation function corresponding to the HN law is $\phi_{\rm HN}(t)=\Pr(\tau_{p}\,S_{\alpha}(\varGamma_{\gamma})^{1/\alpha}\geq t)$. By direct calculations [77] one can easily find

[TABLE]

Using formula (46) for the Laplace transform of completely asymmetric $L\alpha S$ random variables (58), i. e. ${\cal L}(S_{\alpha};t)={\cal L}(Z_{\alpha};t)=e^{-t^{\alpha}}$, together with (62), we come to

[TABLE]

Thus, the effective relaxation rate of the clustered system considered above is

[TABLE]

So, the HN relaxation function can be expressed in the form of a weighted average (31) of an exponential decay $e^{-bt}$ with respect to the distribution $g_{\rm HN}(b)$ of the effective relaxation rate $\widetilde{\beta}_{\rm HN}$. In this case we obtain

[TABLE]

where $\psi(b)=\frac{\pi}{2}-\arctan\left(\frac{(b\tau_{p})^{-\alpha}+\cos(\pi\alpha)}{\sin(\pi\alpha)}\right)$ (Figure 2). At this point we have to stress that such a representation, following the most natural attempt to non-exponential relaxation [41, 10], it is not possible to derive the HN relaxation function with $\gamma>1$. To describe the atypical relaxation data represented in Figure 1, the proposed scheme with random variables $N_{i}$, $M_{j}$, and $\beta_{iN}$ should be modified. Instead of the overestimating the number of clusters and super-clusters (see Eqs.(51) and (52)), the procedure of underestimating these numbers should be involved.

Now the number $K_{N}$ of active entities in the system satisfies the relation

[TABLE]

and the number $L_{N}$ of super-clusters is defined by another dependence

[TABLE]

In this case we analyze (56) in a way analogous to that using the overestimating scheme. The study provided above makes the similar derivations unnecessary, so we omit them and at once write the relaxation rate of the clustered system corresponding to the atypical relaxation data, namely we have

[TABLE]

The subscript JWS has been used to show a link to the relaxation function for atypical relaxation data, derived in the diffusion framework by Jurlewicz, Weron and Stanislavsky (JWS) [78, 79, 80]. The relaxation function, corresponding to the JWS law, has a different form than the HN one [81]. It reads

[TABLE]

Then the pdf of $\tilde{\beta}_{\rm JWS}$ (giving the relaxation function in the form of a weighted average of an exponential decay $e^{-bt}$) takes the very similar (but not the same as $g_{\rm HN}(b)$) form

[TABLE]

where $\psi(b)=\frac{\pi}{2}-\arctan\left(\frac{(b\tau_{p})^{\alpha}+\cos(\pi\alpha)}{\sin(\pi\alpha)}\right)$ (see Figure 2). In relation to Eqs.(5) and (31), for the two-power law dielectric susceptibilities (see Eqs.(3) and (4)), by Tauberian theorem [31], we have the following asymptotic properties of effective relaxation rate pdf

[TABLE]

see the bottom panels in Figure 2. Let us note that the Laplace transform of the generalized arcsine pdf $f_{\gamma}(x)$

[TABLE]

where $M(a,b,x)$ is the Kummer (confluent) function [82], describes a mirror reflection of the CD law in frequency domain as Fig. 2 in [83]. Hence, the random effective waiting time $\widetilde{\theta}_{\rm JWS}$ is given by a mixture of random variables $\tau_{p}\,S_{\alpha}\,(\varXi_{\gamma})^{1/\alpha}$, and the relaxation function takes the form

[TABLE]

This allows us to find the frequency-domain shape function

[TABLE]

where $\varrho_{\gamma}$ is the pdf of the random variable $\Xi_{\gamma}$, i. e.

[TABLE]

for $x>0$. As a result, we get the following expression

[TABLE]

It is useful to compare probabilistic properties of two different random variables, $\varGamma_{\gamma}$ and $\varXi_{\gamma}$, appearing in this context. The distribution of the first of them is characterized by all finite non-zero moments, whereas in the second case all the integer moments become equal to zero. Note that the relaxation patterns very close to the mirror CD law are observed in neo-hexanol ($m=0.72$ and $1-n=0.95$, see Fig. 5.27 in [10]) as well as in gallium (Ga)-doped mixed crystals [84]. The difference among CD, CC and mirror relaxations with the same parameter are shown in Figure 3.

The finiteness of the expected value and long-tail property (11) can be presented only on three different levels: active entity $\Rightarrow$ cluster $\Rightarrow$ cooperative region (super-cluster) of the complex system. To sum up, Table 1 shows the connection between the internal properties of complex system’s dynamics and the empirical relaxation responses, as well as, the physical sense of the parameters characterizing the responses. The proposed approach leads to a very general scenario of relaxation, from the stochastic nature of microscopic dynamics through the hierarchical structure of parallel multi-channel processes to the empirical macroscopic laws of relaxation (see Figure 4).

The fundamental consequence of property (3) is that for large $\omega$ the ratio of the imaginary to real term of the complex susceptibility $\chi(\omega)=\chi^{\prime}(\omega)-i\chi^{\prime\prime}(\omega)$ is a constant, dependent only on the exponent $n$

[TABLE]

The physical significance of this simple property is that at high frequencies the ratio of the macroscopic energy lost per radian to the energy stored at the peak is independent of frequency [10]. However, the D response does not have this property.

A. K. Jonscher [10] has advanced a hypothesis that the fractional power law (3) and the energy criterion (68) are inescapably connected with the fact that the energy loss in every microscopic reversal is independent of the rate of reversals in the corresponding frequency range. He assumed that since in any dielectric system the total polarization is a sum of individual microscopic polarizations and the total loss is the sum of individual microscopic losses, the microscopic relationship also must have the property of energy lost to energy stored being independent of frequency. The fact is based on the identical property of individual structural elements of the systems. This explains the universality in the large scale behavior of complex systems.

Note, the above-mentioned limiting form (57) basically is determined by the tail behavior of $F_{\beta_{i}}(b)$ for large $b$, i. e. by asymptotic properties of $F_{\beta_{i}}(b)$. The detailed knowledge of its other properties is not necessary. The distribution function $F_{\beta_{i}}(b)$ belongs to the domain of attraction of the $L\alpha S$ law with the index of stability $0<\alpha<1$ if and only if [31] for each $x>0$

[TABLE]

This condition can be interpreted as a type of self-similarity:

[TABLE]

The self-similarity was suggested earlier [7, 56, 85] as a fundamental feature of relaxation phenomena. Let us stress that in the limit-theorem approach this result is obtained on the pure probabilistic base, independently of the physical details of systems.

In the framework of the correlated-cluster approach the physical intuition of A. K. Jonscher may be strictly argumentative. Really, the condition (11) applied to any relaxation rate $\beta_{i}$ leads to the scaling property of the relaxation-rate distribution at large $b$ (see also Eq.(70)). The asymptotic behavior of the distribution is connected with the short-time asymptotic properties of the associated relaxation function $\phi(t)$, and the response function as its derivative $f(t)=-d\phi(t)/dt$ takes the form

[TABLE]

for $t\to 0$, where $U(t)$ is a slowly varying function so that $\lim_{t\to 0}U(ct)/U(t)=1$ for any constant $c>0$ [86]. It may be easily verified that the short-time behavior of $f(t)$ corresponds to the high-frequency properties of the susceptibility $\chi(\omega)$:

[TABLE]

The result yields straightforwardly the energy criterion (68) with $n=1-\alpha$. The long-tail property of micro/meso/macroscopic relaxation rates with the parameter $\alpha$ leads to micro/meso/macroscopic energy criterion with the characteristic constant $1-\alpha$. The analysis of the model shows [70, 71, 72] that in the HN, JWS, CC and KWW responses the energy criterion holds for all micro/meso/macro levels, and the power-law exponent $n$ for the HN case (as well as the power-law exponent $m$ for the JWS case) is defined not only by the long-tail property of the distribution of cluster sizes, but also of super-cluster sizes. In the CD case the microscopic energy criterion is not fulfilled. The high-frequency power law of this response results only from the long-tail property of the distribution of super-cluster sizes (see [72] for details).

3.5.3 Finite mean cluster sizes

Assume that the sequences of random variables $N_{i}$, $M_{j}$, and $\beta_{iN}$ are stochastically independent; and each sequence consists of iid non-negative random variables. Assuming moreover finite-average cluster size $\langle N_{i}\rangle=\eta_{0}$ we obtain that $K_{N}\approx\eta_{0}^{-1}N$ (with probability 1) for large $N$. Hence, the random sum $\sum\limits_{j=1}^{L_{N}}M_{j}$ is asymptotically distributed as $\mu_{\gamma}\eta_{0}^{-1}N$, where $\mu_{\gamma}$ is a random variable representing a continuous limit of the random number of randomly sized super-clusters. The random variable $\mu_{\gamma}$ indicates the random space in which the super-clusters exist. If the distribution of $M_{j}$ is heavy-tailed (see Eq.(11)) with the tail exponent $0<\gamma<1$, then the pdf of $\mu_{\gamma}$ is given by $q_{\gamma}(x)=[\Gamma(\gamma)\Gamma(1-\gamma)]^{-1}x^{-1}(x-1)^{-\gamma}$ for $x>1$, and 0 otherwise. Here $\Gamma(\cdot)$ is the gamma function. If the super-cluster average size is finite ($\langle M_{j}\rangle<\infty$), then $\gamma=1$ and $\mu_{1}=1$.

Taking $\beta_{iN}=A_{N}^{-1}\beta_{i}$, where $\beta_{i}$ is independent of the system size $N$ and the normalizing system-size dependent constant $A_{N}$ is the same for each dipole, one can write

[TABLE]

For a finite-average distribution of the individual relaxation rates with $\langle\beta_{i}\rangle=b_{0}$, and for $A_{N}=N$ we obtain

[TABLE]

On the other hand, if the distribution of the individual relaxation rates $\beta_{i}$ is heavy-tailed with the tail exponent $0<\alpha\leq 1$ and the scaling constant equals $b_{0}$, then for the normalizing sequence $A_{N}=(\Gamma(1-\alpha)N)^{1/\alpha}$ we get

[TABLE]

where $S_{\alpha}$ is a completely asymmetric $L\alpha S$ random variable with the index of stability $0<\alpha<1$, independent of $\mu_{\gamma}$; formula (72) coincides with Eq.(71), if we take $\alpha=1$ and $S_{1}=1$.

Now, one can derive the corresponding relaxation or response functions. In particular, the case $0<\gamma<1$ and $\alpha=1$ corresponds to the CD relaxation pattern with $\tau_{p}=\eta_{0}b_{0}^{-1}$, while the case $\gamma=1$ and $0<\alpha<1$ is related to the KWW relaxation function with $\tau_{p}=\eta_{0}^{1/\alpha}b_{0}^{-1}$. Let us remind that $\alpha=1$ refers to $\langle\beta_{i}\rangle=b_{0}<\infty$, and $\gamma=1$ to $\langle M_{j}\rangle<\infty$. For $0<\gamma<1$ and $0<\alpha<1$ one obtains the Generalized Gamma (GG) relaxation for which the relaxation function can be represented as

[TABLE]

where $\varGamma_{\gamma}$ is the gamma distributed random variable with the shape parameter $\gamma$ and the scale parameter equal to 1, and $\tau_{p}=\eta_{0}^{1/\alpha}b_{0}^{-1}$. The corresponding response function $f_{GG}(t)$ takes hence the form of GG pdf [33]

[TABLE]

It should be noticed that this relaxation function is supported in experimental data [87] (propylene glycol and 2-picoline in tri-styrene). As one can see, the GG response results from the heavy-tail properties of both the active-dipole relaxation-rate ($\beta_{iN}$) and the super-cluster-size ($M_{j}$) distributions. Its parameters $0<\alpha,\gamma<1$ are equal to the respective tail exponents. The special case of the CD pattern ($\alpha=1$) corresponds to the finite-average $\beta_{iN}$ distributions taken instead of the heavy-tailed one. Similarly, the KWW response ($\gamma=1$) refers to the finite-average $M_{j}$ distribution. As the generalized model has two spatial scales (clusters and super-cluster regions), one could expect that the latter scale corresponds to larger relaxation times than the scales being related to the clusters. But the clusters of super-cluster regions are dynamically constrained. Therefore, their parameters $\alpha$ and $\gamma$ only influence on short-time behavior of this relaxation via active dipoles. When such a constraint is absent or weak, the clusters and super-cluster regions are responsible for different time scales. Consequently, the super-cluster evolution can determine the relaxation trend in lower frequencies. To sum up, Table 2 demonstrates the connection between the internal properties of such complex system’s dynamics and the parameters characterizing the empirical relaxation responses.

Now we observe that the GG function, as well as KWW and CD ones (see Figure 5), exhibit the short-time power law

[TABLE]

where $n=1-\alpha\gamma$. For the long-time limit both the GG and KWW functions decay stretched exponentially with the exponent $\alpha<1$, while the CD function decays simply exponentially. The time-domain limiting properties of the GG function correspond to those of the frequency-domain response given by

[TABLE]

where ${{}_{1}\!\Psi_{0}}[z\ |\ (\alpha\gamma,\alpha)]=\sum^{\infty}_{n=0}\Gamma(\alpha(\gamma+n))\,z^{n}/n!$ is a special case of the Fox-Wright Psi function [88].

Substituting $s=\omega t$ in (75), we get

[TABLE]

The dielectric susceptibility exhibits hence the fractional high-frequency power law with different fractional exponents for the GG, CD and KWW functions. In the low-frequency range we come to

[TABLE]

while

[TABLE]

Therefore, in the low-frequency limit for the GG pattern (with its special KWW and CD cases) one observes linear dependence on frequency in the absorption term. The results have been obtained in [89]. It should be noticed that the attempt to fit the experimental data analyzed in [87] by two different Havriliak-Negami relaxation functions requires seven various parameters [90] instead of the three ones as in the case [87] based on the generalization of CD and KWW functions. In this context the GG function is more preferable.

3.6 Spatial randomness

The Bernoulli binomial distribution [31]

[TABLE]

is often applied for a finite area study as the exact model with spatial randomness. It is easy to show that when $N\to\infty$ and $p\to 0$ so that $pN$ is a constant, the Bernoulli distribution becomes the Poisson distribution, that describes random patterns in infinitely large areas. Both Bernoulli and Poisson distributions suppose that the probability that an individual of species is found in a given area is independent of the presence of other individuals in the same area. If a space cell already contains an individual, then the cell will be more likely to keep more individuals, whereas empty cells are apt to remain empty. It cannot be considered strictly random from the spatial point of view. The important feature relates to aggregated patterns. To describe such patterns, we use the negative binomial distribution [33]. The latter is often considered as a flexible alternative to the Poisson model for count data. The negative binomial distribution is a substitute for the Poisson distribution, when it is doubtful whether the strict requirements, particularly independence, for the Poisson distribution will be satisfied. In this case the number of entities taking part in the relaxation process is not necessarily fixed with respect to all $N$ dipoles forming the system. The transition process starts with a random initial number $\nu_{M}$ of ordering dipole orientations and then runs due to the individual transitions of dipole orientations at random instants of time $\theta_{1},\theta_{2},\ldots,\theta_{N}$. The number $\nu_{M}$ is an integer-valued random variable depending on the size $N$ of the system. The effective relaxation rate $\widetilde{\beta}$ is given by the summation of individual rates $\beta_{i}$ over all $\nu_{N}$ possible routes for its realization, i. e.

[TABLE]

This form of $\widetilde{\beta}_{N}$, instead of that defined in Eq.(27), yields now a mixed [66] effective relaxation rate $\widetilde{\beta}$. The negative binomial distribution for $\nu_{N}$ is written as

[TABLE]

with $n=1,2,\,\dots$ and the parameter $\gamma>0$. It follows from the book of Johnson and Kotz [33] that in the limit $N\to\infty$ Eq.(78) transforms into the gamma distribution with pdf $g_{\gamma}(y)=y^{\gamma-1}e^{-y}/\Gamma(\gamma)$ for $y>0$. Then the survival probability of the entire system {\rm Pr}\bigl{(}\widetilde{\theta}\geq t\bigr{)}=\left\langle e^{-\widetilde{\beta}t}\right\rangle, where $\widetilde{\beta}$ is the limiting random variable of $\widetilde{\beta}_{N}$ given by (77), can be expressed as a mixture of distributions (see Subsection 2.3)

[TABLE]

what can be identified as the tail of Burr distribution [33, 66]. The parameter $\gamma$ has the same physical (or chemical) sense in both gamma distribution and negative binomial one, i. e. it shows a measure of aggregation in the system. It should be pointed out that the above probabilistic schemes start to work even with $N\approx 10^{5}-10^{6}$. If the aggregation is lack, then

[TABLE]

The limiting case ($\gamma\to\infty$) describes the deterministic number of the $L\alpha S$ contributions (see Eq.(30)). The response function $f(t)$, corresponding to Eq.(79), exhibits the two-power-law asymptotics, namely

[TABLE]

The results are supported in experimental data [91, 92, 93]. Note that by direct calculations the relaxation function (79) leads to the transition rate

[TABLE]

so that for $\gamma\to\infty$ it tends to $\alpha A^{\alpha}t^{\alpha-1}$, the transition rate of KWW relaxation.

Now we discuss interpretation of the obtained model in more details. As it has been established in the framework of the limit theorems of probability theory, the stretched exponential law (30) is the only form of the relaxation decay realized in random distributions (such as Bernoulli and Poisson ones) of species as a null model for the species-area relationship. The hyperbolic law (79) corresponds to spatial aggregations described by the negative binomial distribution. It can arise from a variety of random processes. One of such well-known processes is immigration/birth/death-like scheme. In this case the birth (and death) events are not an independent but a contagious process, meaning that a birth has the tendency to induce more births and a death to induce more deaths. As applied to dielectric relaxation, the birth is a merger of ordered dipole orientations, and the death is their breakup. The macroscopic picture of the dipole evolution is an echo of the spatial distributions. There are parameters standing for both, microscopic and macroscopic dynamics. We have denoted them as $\alpha$ and $\gamma$. The parameter $\alpha$ characterizes a Lévy stable (scaling on different levels) character of random processes participating in the development of such systems. Any complex system consists of many objects. Their statistics is assumed to be a summation procedure of many random variables. The variables and their sum (full or partial) belong to the same probability distribution. This means just its stability. If $\alpha\to 1$, then the transition of the system from an excited state to an equilibrium one looks like nothing else but an exponential decay. Otherwise, when the expected value does not exist, the decay is stretched exponential. The relaxation rates take a continuous value from zero to infinity. Its probability density has a power tail denoting that one cannot avoid very large values of the rates. They make a significant contribution to the relaxation evolution to change its behavior with strongly exponential features into stretched ones. The parameter $\gamma$ detects the dipoles as clustered. The clustering is stochastically independent on the stable character of random processes leading to the stretched decay mentioned above. Therefore, knowing the values of $\alpha$ and $\gamma$, we can answer whether there are clusters or not, as well as clarify the probabilistic character for individual relaxation rates. If $\gamma\to\infty$ and $\alpha\to 1$ (no clusters and no stable distribution of rates), then the relaxation decay is only exponential. In this context the dependence (79) is more general than the stretched exponential output.

3.7 Probabilistic vs deterministic modeling

The concept of time-dependent transition rate $r(t)$ in the study of non-exponential relaxation is not always convenient because of its overloaded form in force of $r(t)=-\frac{d}{dt}\ln\phi(t)$ (see Eqs.(39)- (40), and Eq.(81)), where $\phi(t)$ is the relaxation function under consideration. Moreover, $r(t)$ of the HN case (48) becomes very cumbersome, and it does little in understanding of the relaxation mechanisms. Even in the case of the CC relaxation function the transition rate $r(t)$ is of a special form [94]. It is expressed in terms of a ratio of Mittag-Leffler functions

[TABLE]

where $E_{\alpha,\beta}(x)=\sum_{j=0}^{\infty}x^{j}/\Gamma(\alpha j+\beta)$, and $E_{\alpha}(x)=E_{\alpha,1}(x)$. Although in the kinetic equation (15), the relaxation rate $r(t)$ contains some information about stochastic features of the given relaxing system, this approach is not unique. There exists another view on the kinetic equations. If one gives up the kinetic description based on the derivative of first order and accepts the transition rate as a constant, equal to the material one $\tau_{p}$, then the kinetic equation can be rewritten into another representation (using the fractional operators). For example, the CC relaxation equation takes the following fractional form [95, 96]

[TABLE]

with the initial condition $\phi_{\rm CC}(0)=1$, where $\partial^{\alpha}/\partial t^{\alpha}$ is the Riemann-Liouville fractional derivative [97]. It should be stressed that appearance of the fractional derivative in kinetic equations is caused by the completely asymmetric $L\alpha S$ law, describing the major features of stochastic processes in the relaxation of complex systems. In fact, this suggests the probabilistic interpretation of the fractional calculus. The corresponding pseudo-differential equation for the HN relaxation will be considered below.

4 Relaxation in two-state systems

Relaxation, following the D law, may be simply explained by means of a two-state system [98]. Let $N$ be the common number of entities in a complex system. If $N_{\uparrow}$ is the number of entities in the state $\uparrow$, $N_{\downarrow}$ is the number of entities in the state $\downarrow$ so that $N=N_{\uparrow}+N_{\downarrow}$. Assume that for $t=0$ the system is stated in such an order so that the states $\uparrow$ dominate, namely

[TABLE]

where $n_{\uparrow}$ is the part of entities in the state $\uparrow$, $n_{\downarrow}$ the part in the state $\downarrow$. Denote the transition rates by $B$ defined from microscopic properties of the system (for instance, according to the given Hamiltonian of interaction and the Fermi’s golden rule). In this case the kinetic equation describing the ordinary relaxation (D relaxation) takes the form

[TABLE]

where, as usual, the dotted symbol means the first-order time derivative. The relaxation function for the simple two-state system is

[TABLE]

where $2B=\tau_{p}^{-1}$ is the characteristic material constant. It is easy to see that the steady state of the system corresponds to equilibrium with $n_{\uparrow}(\infty)=n_{\downarrow}(\infty)=1/2$. Clearly its response has an exponential character. However, this happens to be the case when entities relax irrespectively of each other and of their environment. If the entities interact with their environment, and the interaction is complex (or random), their contribution in relaxation already will not result in any exponential decay.

4.1 Stochastic arrow of time

The system of equations (83) is a particular case of the general kinetic equation used for the description of evolution of Markov random processes [57]. It is governed by the deterministic array of time. As it was shown in [99, 100, 101], if one accounts for the subordination of $n_{\uparrow}(\tau)$ and $n_{\downarrow}(\tau)$ by an appropriate random process, the kinetic equation (83) changes its form in dependence of chosen subordinators. Recall, a subordinator itself is a stochastic process (with independent, stationary, non-negative increments) describing the evolution of time within another stochastic process (see Section 2.4). Such a subordinator introduces a new time clock (stochastic time arrow). In fact, the notion of a subordinator allows determining the random number of “time steps”. This concept is very helpful to account for appearance of traps in the evolution of relaxing entities.

Let us consider time evolution of the number of entities in the states $\downarrow$ and $\uparrow$ as parent random processes in the sense of subordination. Assuming that they may be subordinated by another random process with a pdf, say $p(t,\tau)$. Then the entity ratio $m_{\uparrow}(t)$ in the state $\uparrow$ and the entity ratio $m_{\downarrow}(t)$ in the state $\downarrow$ under the temporal subordination are determined by the following integral relations

[TABLE]

To derive the above equation, it is necessary to choose a subordinator in the clear form.

If a new time clock is the inverse $L\alpha S$ subordinator, then the equation of the two-state system (85) takes the following form

[TABLE]

where $D^{\alpha}$ is the $\alpha$-order fractional derivative with respect to time, and $0<\alpha\leq 1$. Here we use the Caputo derivative [102], namely

[TABLE]

where $x^{(n)}(t)$ means the $n$-derivative of $x(t)$. In this case the relaxation function reads

[TABLE]

where $E_{\alpha}(z)$ is the one-parameter Mittag-Leffler function [103]. Evolution of the states $\downarrow$ and $\uparrow$ is shown in Figure 6. Following Eq.(5), the frequency-domain CC function is

[TABLE]

With reference to the theory of subordination [104] the CC law shows that the entities tend to equilibrium via motion alternating with stops so that the temporal intervals between them are random. The random values are just governed by the inverse $\alpha$-stable subordinator.

Notice, by a slight change in the method we can easily obtain not only the CC law of relaxation. The approach has a wide potential that will be presented in the next subsection.

4.2 Memory function

If the distribution of a nonnegative stochastic process $U(\tau)$ is infinitely divisible, then its Laplace transform takes the form

[TABLE]

where $\Psi(s)$ is the Laplace exponent according to the Lévy-Khintchine formula [31] (called also as the Lévy exponent). Many examples of such distributions are well known [33]. Among them are such as Gaussian, inverse Gaussian, $\alpha$-stable, tempered $\alpha$-stable, exponential, gamma, compound Poisson, Pareto, Linnik, Mittag-Leffler and others, including completely asymmetric $\alpha$-stable distributions. Let $f(\tau,t)$ be the pdf of $U(\tau)$. Then the mean $\langle e^{-sU(\tau)}\rangle$ is the Laplace transform of $f(\tau,t)$ equal to

[TABLE]

Knowing $f(\tau,t)$, it is not difficult to find the pdf $g(t,\tau)$ of its inverse process. Based on the Laplace transform of $g(t,\tau)$ with respect to $t$ we come to

[TABLE]

It should be pointed out that the exponential term in the Laplace image $\tilde{g}(s,\tau)$ allows one to simplify further calculations by reducing them to algebraic transformations.

Let us now come back to the two-state system where interactions of entities with their environment are taken into account with help of the temporal subordination under an infinitely divisible processes mentioned above. We will analyze the time evolution of the number of entities in the states $\downarrow$ and $\uparrow$. They are parent random processes in the same manner as in the case of the previous subsection. To account for interaction of entities with the environment, we subordinate the latter processes by another random process with a probability density $g(t,\tau)$ according to the relations (85). In the steady state the system supports the relation $m_{\uparrow}(\infty)=m_{\downarrow}(\infty)=1/2$ conventionally. Then the relaxation function can be written as

[TABLE]

Using the Laplace image $\tilde{g}(s,\tau)$ as applied to infinitely divisible random processes (89), we find the frequency-domain response in a simple form

[TABLE]

where $B$ is a constant transition rate. The latter equation clearly shows that the relaxation behavior will be determined by the Laplace exponent $\Psi(s)$.

In the subordination framework, relaxation of the two-state system is described by the following equations

[TABLE]

where the kernel $M(t)$ plays the role of a memory function. The time-dependent function $M(t)$ can be written as an inverse Laplace transform ${\cal L}_{t}^{-1}$, i. e.

[TABLE]

where $c$ is large enough that $1/\Psi(s)$ is defined for ${\it Re}\ s\geq c$, and $i^{2}=-1$. Then the relaxation function satisfies the following equation

[TABLE]

Using Eq.(92), we get

[TABLE]

The relaxation response $f(t)=-\frac{d\phi(t)}{dt}$ obeys another integral equation

[TABLE]

From the above formula it follows immediately an interesting relationship between $f(t)$ and $M(t)$:

[TABLE]

Note that every infinitely divisible stochastic process is defined by its Laplace (or Lévy) exponent $\Psi(s)$. There exists one-to-one connection between the exponent and the corresponding memory function. Thus, the memory function $M(t)$ is the same important characteristics for relaxing systems as their time-domain relaxation responses. However, it is not easy to detect $M(t)$ in experiments. The physical significance of the memory function consists in its asymptotic properties. Power-law tails of memory functions in short or/and long times reflect directly similar properties in the asymptotic behavior of the corresponding relaxation responses.

4.3 Memory function formalism for empirical laws of relaxation

A significant amount of experimental data on relaxation of the disordered systems supports some types of empirical laws (it is about D, CD, CC and HN) [10, 11]. According to [105], their memory functions take the following forms

[TABLE]

where $0<\alpha,\gamma\leq 1$, $\theta(t)$ denotes the Heaviside step function, and $E_{\alpha,\beta}(x)=\sum_{j=0}^{\infty}\,x^{j}/\Gamma(\alpha j+\beta)$, with $\alpha,\beta>0$, defines the well-known two-parameter Mittag-Leffler function [103, 106]. The memory function connected with the HN law takes a cumbersome form

[TABLE]

expressed in terms of the three-parameter Mittag-Leffler function [75, 107, 108], which has the following Taylor series representation

[TABLE]

The shapes of all these (D, CC, CD and HN) memory functions are shown in Figure 7. They clearly indicate that the memory function $M_{\rm HN}(t)$ is the most general case leading to particular cases $M_{\rm D\ }(t)$, $M_{\rm CC}(t)$ and $M_{\rm CD}(t)$ in limit values of parameters ($\alpha\to$ 1 and/or $\gamma\to$ 1). Really, based on relation $E^{\mu}_{1,\mu}(y)=e^{y}/\Gamma(\mu)$, the HN memory function easily transforms into the CD ($\alpha=1$) and D ($\alpha=\gamma=1$) memory ones. When $\gamma=1$, we obtain the CC memory function using the sum $\sum^{\infty}_{j=0}\tau_{p}^{-\alpha j}\,t^{\alpha j}\,E^{j+1}_{\alpha,\alpha(j+1)}(-\tau_{p}^{-\alpha}t^{\alpha})=1/\Gamma(\alpha)$.

5 Anomalous diffusion approach. Stochastic processes

Any relaxation process is seen as a change in time (growth or decay) of a macroscopic physical magnitude characteristic for the observed system (e.g. polarization-depolarization of a dielectric material) is accompanied by diffusion of a corresponding physical variable. In particular, the decay or growth of polarization is accompanied by diffusion of dipole orientations termed as diffusion of an excitation mode. From the theoretical point of view, the models should be adequate to the dual description of the kinetic process leading however to the same results, i.e. to the empirical evidence. Taking into account the stochastic nature of physical mechanisms underlying the observed time evolution of the system, the models from one side should yield a rigorous definition of the probability of changing or not the initial state of the system and, from the other side, the models should yield a rigorous stochastic process representing the underlying diffusion. Only this dual description can give full information on anomalous dynamical properties of various systems.

In this context the non-exponential relaxation of complex systems can be modeled from the idea based on an excitation undergoing diffusion in the systems [8, 9]. Then the relaxation function $\phi(t)$ is determined by the temporal decay of a given mode $k$, and in the framework of the one-dimensional continuous time random walks (CTRWs) it is given by the following Fourier transform (supported on whole the real axis)

[TABLE]

where $k>0$ has the physical sense of a wave number, and $\widetilde{R}(t)$ denotes the diffusion front (scaling limit of the CTRW) under consideration. If $\widetilde{R}(t)$ is supported on the positive half-line, its Fourier transform is replaced by the Laplace transform [109]. Thus we obtain

[TABLE]

The details will be given in the next subsections.

5.1 Uncoupled continuous time random walk. Diffusion. Cole-Cole relaxation

The idea of the continuous time random walk (CTRW) was firstly proposed by E. W. Montroll and G. H. Weiss in 1965 [110]. Although the term “random walk” was introduced by K. Pearson in 1905 [111], the formalism of simple random walks was known else in the XVII-th century. The random walk approach is based on assumption that step changes are made through equal time intervals. This was the first (and rough) approximation for modeling many various physical, chemical and economical phenomena [112]. The CTRW model went further to study a random waiting time among subsequent random jumps.

Briefly, the core of the CTRW methodology is the following. Consider a sequence ${T_{i}}$, $i=1,2,\dots$ of non-negative iid random variables which represent waiting-time intervals between subsequent random jumps $R_{i}$ of a walker. The random time interval of $n$ jumps equals

[TABLE]

and the space position of the walker after the $n$ jumps is given by the sum

[TABLE]

where $R_{i}$ are real iid variables showing both the length and the direction of the i-th jump. The random variables $R_{1},R_{2},\dots$ are assumed to be independent of $T_{1},T_{2},\dots$, although this is not obligatory. The variables $R_{i}$ may be multi-dimensional vectors. Without loss of generality, we consider the walks to be one-dimensional.

The random number $N_{t}$ of jumps, performed by a walker up to time $t>0$, can be determined by the largest index $n$ for which the sum of $n$ interjump time intervals does not exceed the observation time $t$, namely

[TABLE]

The process $N_{t}$ is often called the counting process or otherwise as the renewal process. The total distance attaining by the walker after the $N_{t}$ jumps becomes then

[TABLE]

The cumulative stochastic process (100) is just known as the CTRW.

Though the aforesaid walks consist only of discrete time and space steps, the random walk model can be generalized to “continuous steps” being, hence, closer to the physical reality. The label continuous indicates the fact that the index $t$ in the CTRW belongs to a continuous set $[0,\infty)$, but does not imply the continuity of the paths. Assuming that the interjump time intervals $T_{i}$ belong to the domain of attraction of a completely asymmetric $L\alpha S$ distribution with the index $0<\alpha<1$ ($\beta=1$), i. e.

[TABLE]

Then the limit theorems [113, 114] yield the continuous limit of the random sum (97), i. e.

[TABLE]

where $c^{-1/\alpha}$ is the time-rescaling constant chosen appropriately, $\lfloor x\rfloor$ denotes the integer part of $x$, and $U_{\alpha}(\tau)$ is a strictly increasing $L\alpha S$ process [104]. Let the jumps $R_{i}$ belong to the domain of attraction of a $L\eta S$ distribution $S_{\eta,\beta}(x),0<\eta\leq 2,|\beta|\leq 1$ so that the continuous limit reads

[TABLE]

where $X_{\eta}(\tau)$ is a $L\eta S$ process known as the parent process. If $\eta=2$ ($\beta=0$), the parent process $X_{2}(\tau)$ becomes the classical Brownian motion. Both, the process $U_{\alpha}(\tau)$ and the process $X_{\eta}(\tau)$ are indexed by random operational (internal) time $\tau$.

As $T(n)$ is the random time interval elapsed among $n$ jumps, and $N_{t}$ is the number of jumps occurred up to time $t>0$, they are connected with each other by the following formula

[TABLE]

Using Eq.(102), we get

[TABLE]

where $S_{\alpha}(t)$ is the inverse $L\alpha S$ subordinator, relating the internal and the observable times, and $c\to\infty$. In fact, the continuous limit of the discrete counting process $N_{t}$ is the hitting time process $S_{\alpha}(t)=\inf\{\tau\!:\!U_{\alpha}(\tau)>t\}$ satisfying the relation $S_{\alpha}(U_{\alpha}(\tau))=\tau$ (almost surely and almost everywhere). Therefore, this process can be treated as the inverse to $U_{\alpha}(\tau)$. Because the process $U_{\alpha}(\tau)$ is strictly increasing, the process $S_{\alpha}(t)$ is nondecreasing. The hitting time $S_{\alpha}(t)$ is called also a first passage time [31]. For a fixed time it represents the first passage of the stochastic time evolution above this time level. Note, the random process $S_{\alpha}(t)$ just depends on the true time $t$, and any sample trajectory of the process $S_{\alpha}(t)$ can be only increasing. Both processes $X_{\eta}(\tau)$ and $U_{\alpha}(\tau)$ depend on $\tau$. To find the position of a walker at the real (true) time in this irregular motion, we derive the continuous limit of Eq.(100), i. e.

[TABLE]

known as the anomalous diffusion process (diffusion front $\widetilde{R}(t)$) according to, for example, [115]), directed by the inverse $L\alpha S$ subordinator $S_{\alpha}(t)$. In this issue the process $S_{\alpha}(t)$ plays the role of a new time clock (stochastic time arrow) in $X_{\eta}(\tau)$ instead of $\tau$. In the framework of this approach the stochastic process (105) opens perspectives to describe an arbitrary diffusion.

Let us discuss a relationship between the pdf of the position $r_{t}$ of a walking particle at real time $t$ and the probability densities of random processes $X_{\eta}(\tau)$ and $S_{\alpha}(t)$. As applied to $X_{\eta}(S_{\alpha}(t))$ defined in Eq.(105), we can consider it as a subordination of processes. If the processes $X_{\eta}(\tau)$ and $U_{\alpha}(\tau)$ are uncoupled (i. e. independent on each other), the pdf of $r_{t}$ with $t\geq 0$ can be written as a weighted integration over the internal time $\tau$ so that

[TABLE]

where $f(x,\tau)$ denotes the pdf to find the parent process $X_{\eta}(\tau)$ at $x$ on operational time $\tau$ and $g(t,\tau)$ is the pdf to be at the operational time $\tau$ on real time $t$. Then the Fourier transform $\widehat{p}(k,t)=\langle\exp(ikX_{\eta}[S_{\alpha}(t)])\rangle$ and the Laplace transform $\widetilde{p}(k,t)=\langle\exp(-kX_{\eta}[S_{\alpha}(t)])\rangle$ take the following forms

[TABLE]

where $k>0$ is the wave number mentioned above. Consequently this allows one to determine the relaxation function from the probabilistic properties of diffusive processes (see Eqs.(95) and (96)).

The Laplace image for pdf of a non-negative $L\alpha S$ variable is

[TABLE]

where $0<\alpha<1$. If $f(\tau,t)$ is the pdf of $U_{\alpha}(\tau)$, then the pdf $g(t,\tau)$ of its inverse $S_{\alpha}(t)$ reads

[TABLE]

Taking the Laplace transform of $g(t,\tau)$ with respect to $t$, we get

[TABLE]

It follows from Eq. (110) that the pdf of the inverse $L\alpha S$ process is

[TABLE]

where $Br$ denotes the Bromwich path [82], and the function $F_{\alpha}(z)$ is given by the following Taylor series expansion

[TABLE]

For more details see [95, 101]. In this framework, the relaxation function $\phi(t)$ that describes the temporal decay of a given mode $k$, can be expressed through the Fourier transform (95) of the diffusion process $X_{\eta}(S_{\alpha}(t))$. Consequently, this case leads to the CC relaxation function $\phi_{CC}(t)=E_{\alpha}\left[-(t/\tau_{p})^{\alpha}\right]$, where $\tau_{p}\sim|k|^{-\eta/\alpha}$ [100, 104]. The frequency-domain description of the CC relaxation takes the form (87). If $\alpha=1$, it becomes the classical D relaxation.

5.2 Coupling between the very large jumps in physical and operational times

Let us now make use of the fact that the stochastic time evolution $U_{\alpha}(\tau)$ and its (left) inverse process $S_{\alpha}(t)$ permits one to underestimate or overestimate the time interval $T(N_{t})$ of $N_{t}$ random steps performed up to the physical time $t$ at which the position of a walker is observed:

[TABLE]

what follows directly from the definition (97). In fact, the two processes $T(N_{t})$ and $T(N_{t}+1)$ correspond to underestimating and overestimating the real time $t$ from the random time steps $T_{i}$ of the CTRWs.

In terminology of the Feller’s book [31] the variable

[TABLE]

is the residual waiting time (life-time) at the epoch $t$, and

[TABLE]

is the spent waiting time (age of the entity that is alive at time $t$). Importance of these random variables can be explained by one remarkable property. For $t\to\infty$ the variables $Y_{t}$ and $Z_{t}$ have a common, proper limiting distribution only if their probability distributions $F(y)$ and $F(z)$ have finite expectations. However, if the distribution $F(x)$ satisfies

[TABLE]

where $0<\alpha<1$ and $L(xt)/L(t)\to 1$ as $x\to\infty$, then according to [116], the pdf of the normalized variable $Y_{t}/t$ is given by the generalized arcsine law

[TABLE]

while $Z_{t}/t$ obeys

[TABLE]

Since $\Sigma_{N_{t}}=t-Y_{t}$ and $\Sigma_{N_{t}+1}=Z_{t}+t$, the distributions of $\Sigma_{N_{t}}/t$ and $\Sigma_{N_{t}+1}/t$ can be obtained from Eqs. (112) and (113) by a simple change of variables $1-x=y$ and $1+x=z$, respectively.

In this case $T(N_{t})/t$ tends in distribution in the long-time limit to random variable $Y$ with density

[TABLE]

and $T(N_{t}+1)/t\stackrel{{\scriptstyle d}}{{\rightarrow}}Z$ with the pdf equal to

[TABLE]

Both pdfs, $p^{Y}(x)$ and $p^{Z}(x)$, are special cases of the well-known beta density. It should be noticed that the density $p^{Y}(x)$ concentrates near 0 and 1, whereas $p^{Z}(x)$ does near 1. Near 1 both densities tend to infinity. This means that in the long-time limit the most probable values for $T(N_{t})$ occur near 0 and 1, while for $T(N_{t}+1)$ they tend to be situated near 1.

The nonequality (111) can also be represented in a schematic picture of time steps. Next, a passage from the discrete process $T_{i}$ to the continuous one $U_{\alpha}(\tau)$ allows one to reformulate the unequality (111) as

[TABLE]

underestimating or overestimating the real time $t$. Note, the difference $U_{\alpha}(S_{\alpha}(t))-t$ is the leap-over processes [117, 118]. The pdfs of $U^{-}_{\alpha}(S_{\alpha}(t))$ and $U_{\alpha}(S_{\alpha}(t))$ take the forms, respectively,

[TABLE]

valid for any time $t>0$ (see Figure 8). The moments of $U^{-}_{\alpha}(S_{\alpha}(t))$ and $U_{\alpha}(S_{\alpha}(t))$ can be calculated directly from the moments of $Y$ and $Z$ by using the following relations

[TABLE]

Thus, the process $U^{-}_{\alpha}(S_{\alpha}(t))$ has finite moments of any order, while $U_{\alpha}(S_{\alpha}(t))$ gives us even no finite the first moment. The process $U_{\alpha}(S_{\alpha}(t))>t$ is too long also in the limit formulation. Notice that $p^{+}(t,y)=y^{-2}p^{-}(t^{-1},y^{-1})$. In our construction of the compound subordinators $U^{-}_{\alpha}(S_{\alpha}(t))$ and $U_{\alpha}(S_{\alpha}(t))$ the processes $U_{\alpha}(\tau)$ and $S_{\alpha}(t)$ are clearly coupled.

5.3 Anomalous diffusion with under- and overshooting subordination

According to studies presented in [79], the widely observed fractional two-power relaxation dependencies (3) and (4) are closely connected with the under- and overshooting subordination

[TABLE]

where

[TABLE]

the processes $Y^{U}_{\gamma}(t)$ and $Y^{O}_{\gamma}(t)$ being nothing else as $U^{-}_{\gamma}(S_{\gamma}(t))$ and $U_{\gamma}(S_{\gamma}(t))$ with the index $\gamma$ (Figure 9a). Their main feature is that they are subordinated by an independent inverse $L\alpha S$ process $S_{\alpha}(t)$ forming the compound subordinators $Z^{U}_{\alpha,\gamma}(t)$ and $Z^{O}_{\alpha,\gamma}(t)$, respectively (Figure 9b). The approach enlarges the class of diffusive scenarios in the framework of the CTRWs. Assuming a heavy-tailed cluster-size distribution with the tail exponent $0<\gamma<1$, the coupling between jumps and interjump times tends to the compound operational times $Z^{U}_{\alpha,\gamma}(t)$ and $Z^{O}_{\alpha,\gamma}(t)$ as under- and overshooting subordinators, respectively.

The overshooting subordinator yields the anomalous diffusion scenario leading to the well-known HN relaxation pattern (1), and the undershooting subordinator leads to a new relaxation law (67). Derivation of both cases is presented below (see Section 5.4). These results are in agreement with the idea of a superposition of the classical (exponential) Debye relaxations [41]. Let $X_{\eta}(\tau)$ with $0<\eta\leq 2$ be the parent process that is subordinated either by $Z^{U}_{\alpha,\gamma}(t)$ or $Z^{O}_{\alpha,\gamma}(t)$. Then the subordination relation, expressed by means of a mixture of pdfs, takes the form

[TABLE]

where $p^{\,r}(x,t)$ is the pdf of the subordinated process $X_{\eta}[Z^{U}_{\alpha,\gamma}(t)]$ (or $X_{\eta}[Z^{O}_{\alpha,\gamma}(t)]$) with respect to the coordinate $x$ and time $t$, $p^{X}(x,\tau)$ the pdf of the parent process, $p^{\pm}(y,\tau)$ the pdf of $U^{-}_{\gamma}(S_{\gamma}(t))$ and $U_{\gamma}(S_{\gamma}(t))$ respectively, and $p^{S}(t,\tau)$ the pdf of $S_{\alpha}(t)$. Note, the processes $X_{\eta}[Z^{U}_{\alpha,\gamma}(t)]$ and $X_{\eta}[Z^{O}_{\alpha,\gamma}(t)]$ can be factorized in such a way, that if $X_{2}(\tau)$ is the Brownian process, then, e. g. (see Eq.(105)), the diffusion front $X_{2}(S_{\alpha}(t))$ becomes a mixture of the Gaussian and completely asymmetric $L\alpha S$ laws. This issue will be considered in more details in Section 7 (see also [109, 79]).

The subordinator $U_{\gamma}(S_{\gamma}(t))$ results in stretching of the real time $t$. It will underline scaling properties in short and long times, respectively. Namely, the interesting feature is observed in both CD and HN relaxations. Let us consider the process $U_{\gamma}(S_{\gamma}(t))$ as a subordinator indexed by $\gamma$, i. e. the process is obtained from a $L\gamma S$ random process. Following the temporal decay (95) of a given mode $k$, the relaxation function takes the form

[TABLE]

Here the subscript CD is not by chance. It shows a direct connection of the relaxation function with the CD law [10]. In fact, the one-sided Fourier transform (5) gives

[TABLE]

It should be mentioned that the method of subordination suggests also one more scenario leading to the CD relaxation. It is based on the inverse tempered $L\alpha S$ process (see [119] in more details) that will discuss below.

5.4 Relaxation process from compound subordinators

As it has been already shown above, the CC and CD relaxations are only special cases of the more general HN law. To get that law, the operational time $S_{\alpha}(t)$ of the CC diffusion mechanism has to be modified [79] by means of coupling between jumps and interjump times in the underlying decay of a given mode, representing excitation undergoing diffusion in the relaxing system. The relaxation response will be characterized by short- and long-time power laws with different fractional exponents (as in the HN case) only if the anomalous diffusion scenario is based on a compound operational time. To construct such an operational time, denote conveniently the discussed above processes $U^{-}_{\gamma}(S_{\gamma}(t))$ and $U_{\gamma}(S_{\gamma}(t))$ as $Y^{U}_{\gamma}$ and $Y^{O}_{\gamma}$, respectively. They corresponds to the under- and overshooting subordination scenarios [79]. Next, we can write $Z^{U}_{\alpha,\gamma}(t)\leq S_{\alpha}(t)\leq Z^{O}_{\alpha,\gamma}(t)$ for $t\geq 0$, where $Z^{U}_{\alpha,\gamma}(t)=Y^{U}_{\gamma}[S_{\alpha}(t)]$, $Z^{O}_{\alpha,\gamma}(t)=Y^{O}_{\gamma}[S_{\alpha}(t)]$. The overshooting subordinator $Z^{O}_{\alpha,\gamma}(t)$ leads (stretching the operational time $S_{\alpha}(t)$) to the HN relaxation in the form

[TABLE]

By direct calculations of the Fourier transformation (5) the frequency-domain shape function reads

[TABLE]

The above approach demonstrates clearly a success in probabilistic treatment of the observed relaxation laws (see Figure 10). Therefore, we continue our analysis as applied to the undershooting (compressing $S_{\alpha}(t)$) subordinator $Y^{U}_{\gamma}[S_{\alpha}(t)]$.

In this case we obtain the JWS relaxation function

[TABLE]

which can be identified as a special case of the three-parameter Mittag-Leffler function [107, 108]

[TABLE]

where $(\gamma,j)=\gamma(\gamma+1)(\gamma+2)\dots(\gamma+j-1)$ is the Appell’s symbol with $(\gamma,0)=1$, $\gamma\neq 0$. To avoid any confusion, it should be mentioned that the two-parameter Mittag-Leffler function $E_{\alpha,\beta}(x)$, more common in literature [103], is a special case of the function $E_{\alpha,\beta}^{\gamma}(x)$ with $\gamma=1$. From the series expansion of the simplest Mittag-Lefller function $E_{\alpha}(x)$ it is easy to check by direct calculations of Eq.(122) that

[TABLE]

This type of the relaxation function as an integral has been derived in the CTRW framework [78, 79], and the exact functional form has been established in [80]. Using the relationship with the Mittag-Leffler function, we may now write the kinetic equation for (122) in the pseudodifferential equation form

[TABLE]

where $\partial^{\alpha}/\partial t^{\alpha}$ is the Riemann-Liouville fractional derivative [120], and $\phi_{\rm JWS}(0)=1$ the initial condition. Taking the Fourier transform (5), we get the shape function corresponding to (122) in the form useful for fitting the atypical dielectric spectroscopy data

[TABLE]

If one multiplies both numerator and denominator of the fraction by $(i\omega/\omega_{p})^{\alpha\gamma}$, this clearly leads to Eq.(67). This result points to the following relationship with the HN relaxation function

[TABLE]

For $\gamma\to 1$ the pdf of $Y^{U}_{\gamma}$ tends to the Dirac $\delta$-function. It is easy to check $\lim_{\gamma\to 1}\phi_{\rm JWS}(t)=\phi_{\rm CC}(t)$, and the kinetic equation for $\phi_{\rm JWS}(t)$ takes the mentioned above CC form (82).

It is interesting to compare the power-law characteristics of the relaxation response for both overshooting and undershooting schemes of anomalous diffusion. For both scheme the exponents $n$ and $m$ fall in the range $(0,1)$, and are defined by the subordinator parameters $0<\alpha<1$ and $0<\gamma<1$. The HN relaxation, resulting from the undershooting scheme, is characterized by the exponents $m=\alpha$ and $m>1-n=\alpha\gamma$, and it fits the so-called typical dielectric spectroscopy data, while the JWS case, resulting from the overshooting scheme, demonstrates up-down with $m=\alpha\gamma$ and $m<1-n=\alpha$ and it fits the atypical two-power-law relaxation pattern which, as shown by the experimental evidence (see Figure 1), cannot be neglected (see e. g. [10, 11, 30] and references therein). Such an atypical behavior has been also observed by us in gallium (Ga)-doped Cd0.99Mn0.01Te mixed crystals [121], where sample frequency-domain data measured for Cd0.99Mn0.01Te:Ga at 77K is fitted with the function (67). This material belongs to the semiconductor of group II-VI possessing deep metastable recombination centers. Formation of such centers in Cd1-xMnxTe:Ga results from the bistability of Ga dopant which makes this mixed crystal as an attractive material for holography and high-density data storage (optical memories).

In the language of subordinators we may conclude that the process $Y^{U}_{\gamma}(t)$ makes a rescaling for small times, and the process $Y^{O}_{\gamma}(t)$ turns on a similar rescaling for long times. As for $0<\alpha,\gamma\leq 1$, in the case of HN relaxation the declination of the imaginary susceptibility $\chi^{\prime\prime}(\omega)$ for low frequencies will be greater than for high frequencies, whereas the atypical relaxation shows an opposite relation (see Figure 2).

The original HN relaxation [10, 11] with exponents $0<\alpha,\gamma\leq 1$ satisfies $m\geq 1-n$. Its modified version [30], proposed to fit relaxation data with power-law exponents satisfying $m<1-n$, assumes $0<\alpha,\alpha\gamma\leq 1$. Unfortunately, the HN function with $\gamma>1$ cannot be derived within the framework of diffusive relaxation mechanisms. Only for $\gamma\leq 1$ the origins of the HN function can be found within the fractional Fokker-Planck [122] and the CTRW [79] approaches. The approach considered above includes all the data in the mathematically unified approach.

5.5 Clustered continuous time random walk. Compound counting process

The CTRW process $R(t)$ determines the total distance reached by a random walker until time $t$. It is characterized by a sequence of iid spatio-temporal random steps $(R_{i},T_{i}),i\geq 1$. If we assume stochastic independence between jumps $R_{i}$ and waiting times $T_{i}$, we get a decoupled random walk; otherwise we deal with a coupled CTRW. The distance reached by the walker at time $t$ is given by the following sum

[TABLE]

where $N_{t}=\max\{n:\,\sum_{i=1}^{n}T_{i}\leq t\}$, counts the steps performed up to $t>0$.

Theoretical studies of the relaxation phenomenon in the above framework are based on the idea of an excitation undergoing (anomalous, in general) diffusion in the system under consideration [123]. The relaxation function $\phi(t)$ is then defined by the inverse Fourier or the Laplace transform (see Eqs.(95) and (96)) of the diffusion front $\widetilde{R}(t)$ which is

[TABLE]

where the dimensionless rescaling parameter $t_{0}\to 0$ and $f(t_{0})$ is appropriately chosen renormalization function. The diffusion front $\widetilde{R}(t)$ approximates a position at time $t$ of the walker performing rescaled spatio-temporal steps $(R_{i}/f(t_{0}),t_{0}T_{i})$. The characteristics of the relaxation process are related to the properties of the diffusion front resulting from assumptions imposed on the spatio-temporal steps of the random walk. For example, the decoupled CTRW with power-law waiting-time distributions (i. e., with the random variables $T_{i}$ satisfying $\Pr(T_{i}\geq t)\sim(t/\tau_{0})^{-\alpha}$ as $t\to\infty$ with some $0<\alpha<1$ and $\tau_{0}>0$) leads to the CC relaxation (see Subsection 5.1). But the frequency-domain CC relaxation with the corresponding time-domain Mittag-Leffler pattern is only one of the cases measured in various experiments with complex media, and derivation of those more general patterns requires considering diffusion scenarios based on a compound coupled CTRW representation. It should be pointed out that the simple coupling of type $R_{i}\sim T^{p}_{i}$ (with positive power exponent $p$) does not lead, however, behind the CC relaxation [109, 124]. In contrast, introducing a dependence between the jumps and waiting times by a random clustering procedure (being a stochastic generation of the well-known deterministic renormalization approaches) we can derive another relaxation laws like the CD or HN ones [78, 79, 80]. Below, we present the clustered CTRW scenario which leads directly to the results discussed in the preceding section.

Let $M_{j}$ be a sequence of iid positive integer-valued random variables independent of the pairs $(R_{i},T_{i})$. Next, assume that the jumps and waiting times are assembled into clusters of random sizes $M_{1},M_{2},\dots$. This assumption allows one to transform the sequence of spatio-temporal steps $(R_{i},T_{i})$ into a new sequence $(\overline{R}_{j},\overline{T}_{j})$ of random sums

[TABLE]

Then the position $R^{M}(t)$ of the walker is determined by $(\overline{R}_{j},\overline{T}_{j})$ and, in accordance with the general formula (123), it is given by

[TABLE]

where $\bar{\nu}(t)=\max\{n:\,\sum_{j=1}^{n}\overline{T}_{j}\leq t\}$. The dependence between the jumps $\overline{R}_{j}$ and the waiting times $\overline{T}_{j}$ of the coupled CTRW process $R^{M}(t)$ is determined by the distribution of the cluster sizes $M_{j}$.

As an example, let us consider the simple case of a random walk, when the waiting times are represented by equal intervals in time, i. e. $T_{i}=\Delta t$. In this case we have

[TABLE]

where the step-clustering procedure (124) yields then $\overline{T}_{j}=M_{j}\Delta t$ for $j\geq 1$, and the coupled process $R^{M}(t)$ in (125) takes the following form

[TABLE]

Here $U^{M}(\bar{\nu}(t))$ is a compound counting process obtained from

[TABLE]

and

[TABLE]

Observe that formula (127) is an analog of (126) with the compound counting process $U^{M}(\bar{\nu}(t))$ substituting the deterministic number $\lfloor t/\Delta t\rfloor$ of performed jumps $R_{i}$. The counting process $U^{M}(\bar{\nu}(t))$ is always less than $\lfloor t/\Delta t\rfloor$, and it is hence a special case of the undershooting compound counting process [79]. It is also a clear signature of the spatio-temporal coupling provided by the clustering procedure (124).

The idea of compound counting processes in CTRW approach is not new in physics. The resulting CTRW processes were examined in the context of the rareness hypothesis in the fractal-time random walk models (see, e. g. [125, 126]). In general, the compound counting process cumulates random number of random events. Physical situations where the relevance of this scheme holds are numerous. For instance, they take into account the random magnitude of claims’ sequence in insurance risk theory [127], the energy release of individual earthquakes in geophysics or random water inputs flowing into a dam in hydrology [128] where summing the individual contributions yields the total amount of the studied physical magnitude over certain time intervals.

In the classical waiting-jump CTRW idea [110], in which the jump $R_{i}$ occurs after the waiting-time $T_{i}$ (Figure 11a), the random number of the particle jumps performed by time $t>0$ is given by the renewal counting process (99). Then the location of a particle at time $t$ is given by the random sum

[TABLE]

In the alternative jump-waiting CTRW scenario (Figure 11b), the particle jump $R_{i}$ precedes the waiting time $T_{i}$. 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

[TABLE]

This is called the overshooting CTRW, or briefly OCTRW, see [129]. In summary, the CTRW process $R^{-}(t)$ and the OCTRW process $R^{+}(t)$ are obtained by 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, the subordination modifies a random process, replacing the deterministic time index by a random clock process, which usually represents a second source of uncertainty. 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 types of relaxation, are the same. On the other hand, if the coupled case (i. e., dependent coordinates in the random vector $(\overline{T}_{i},\overline{R}_{i})$) is considered, the waiting-jump and jump-waiting schemes may lead to essentially different relaxation patterns.

An important and useful example of the coupled CTRW, different from the most popular Lévy walk [130], has been identified using the clustered CTRW concept, introduced in [131] and developed further in [132]. While in the Lévy walk the jump size is fully determined by the waiting time $R_{i}\sim T^{p}_{i}$ (or equivalently, by flight duration), in the clustered CTRW coupling arises from random renormalization of the number of jumps. The clustered CTRW scheme is relevant to numerous physical situations, including the mentioned above 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 [125, 126, 133, 134]. In all these cases, summation of the individual contributions yields the total amount (in general, random) of the studied physical magnitude over certain time intervals. In the clustered CTRW, the waiting time and the subsequent jump are both random sums with the same random number of summands, and as a consequence, this type of the CTRW is coupled, even if the original CTRW before clustering had no dependence between the corresponding waiting times and jumps (Figure 12). 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 limits, and their governing equations [131, 132].

Considering those cases crucial for modeling of the relaxation phenomena, assume that the waiting times $T_{i}$ have a heavy-tailed distribution with parameter $0<\alpha<1$; i. e., for some $\tau_{0}>0$ we have

[TABLE]

for large $t$. Moreover, let the jump distribution be symmetric and belong to the normal domain of attraction of a symmetric $L\eta S$ law with the index of stability $\eta$, $0<\eta\leq 2$; i. e., for some $\rho_{0}>0$ we have

[TABLE]

for large $x$ if $0<\eta<2$ or $0<\langle{\rm D}^{2}R_{i}\rangle=\rho_{0}^{2}<\infty$ if $\eta=2$. Then, a heavy-tailed distribution of cluster sizes with tail exponent $0<\gamma<1$ yields different anomalous diffusion limits

[TABLE]

and

[TABLE]

of clustered CTRW and clustered OCTRW, respectively [79] and [132]. Here $C$ is a positive constant dependent on the tail exponents $\alpha$, $\eta$ and proportional to the scaling parameter $\rho_{0}$. The parent process $X_{\eta}(\tau)$ is a symmetric $L\eta S$ process (recall, for $\eta=2$ it is just the standard Brownian motion). The directing process $S_{\alpha}(t)$ is defined as

[TABLE]

where $U_{\alpha}(t)$ is a strictly increasing $L\alpha S$ process with the stability index $\alpha$. The undershooting $Y_{\gamma}^{U}(t)$ and overshooting $Y_{\gamma}^{O}(t)$ processes have the same distributions as the processes $U_{\alpha}^{-}(S_{\alpha}(t))$ and $U_{\alpha}(S_{\alpha}(t))$, respectively, for the value of parameter $\alpha$ equal to $\gamma$.

For the clustered CTRW and OCTRW models, the anomalous diffusion fronts $\widetilde{R}^{\pm}(t)$ in (130) and (131) have frequency-domain shape functions that can be identified respectively as the JWS and HN functions, see [79, 80]. The characteristic material constant $\omega_{p}$, appearing in both functions, takes the form

[TABLE]

The diffusion fronts (130) and (131) allow us to develop analytical formulas (see Eq.(95)) for the corresponding time-domain relaxation functions $\phi_{\rm JWS}(t)$ and $\phi_{\rm HN}(t)$ in terms of the three-parameter Mittag-Leffler function [107, 108].

Observe that

[TABLE]

where $p^{\pm}(x,t)$ are pdfs of the anomalous diffusion processes $\widetilde{R}^{\pm}(t)$. Hence, the Fourier-Laplace (FL) images ${\cal I}_{\bf FL}(p^{\pm})(k,s)$ of the functions $p^{\pm}(x,t)$ are just the Laplace images of the corresponding time-domain JWS and HN relaxation functions. According to the results in [132, 135, 136], we can derive $p^{-}(x,t)$ as a mild solution (see about mild solutions of differential equations in [137] for more details) of the following fractional pseudo-differential equation

[TABLE]

where $\delta(x)$ is the Dirac delta function. Equivalently,

[TABLE]

and hence

[TABLE]

with $\omega_{p}^{\alpha}=|Ck|^{\eta}/\tau_{0}^{\alpha}$. Using the following Laplace transformation [107]

[TABLE]

we obtain

[TABLE]

Similarly, the pdf $p^{+}(x,t)$ is a mild solution of equation

[TABLE]

where $p^{X}(x,t)$ and $p^{S}(t,\tau)$ are the pdfs of $X_{\eta}(t)$ and $S_{\alpha}(t)$, respectively [132, 135]. This implies that

[TABLE]

As a consequence, we have

[TABLE]

and

[TABLE]

The short- and long-time behaviors of these functions, which exactly follow the high- and low-frequency power laws (see Subsection 1.2), read

[TABLE]

and

[TABLE]

Thus, the clustered CTRWs and their diffusion limits illuminate the role of random processes in the parametrization of relaxation phenomena. The $L\eta S$ parent process $X_{\eta}(\tau)$, that models particle jumps, does not change the exponents of the relaxation power laws. It only affects the material constant $\tau_{p}$ by determining the spatial features of the anomalous diffusion $\widetilde{R}^{\pm}(t)$. The index $\alpha$ of the process $S_{\alpha}(t)$ that codes the waiting times between jumps, and the index $\gamma$ of the clustering process $Y_{\gamma}(t)$, determine the power-law behavior of the relaxation function in time and frequency. These coefficients characterize the complex dynamics of relaxing systems.

6 Extensions

In real situations there are sufficiently many factors truncating the distribution of $L\alpha S$ processes [138, 139, 140, 141, 142]. This leads to the tempered $L\alpha S$ processes which have all the moments. As it is shown in [119], from the subordination by the inverse tempered $L\alpha S$ process one can derive the tempered diffusion equation and the relaxation function describing the D, CC and CD types of relaxation. The tempered diffusion has a transient character, i. e. a crossover from subdiffusion at short times to normal diffusion at long times. The transient subdiffusion has impact on kinetics of magnetic bright points on the Sun [143] and has been observed in cells and cell membranes [144, 145, 146]. Physical arguments for appearance of such effects are that subdiffusion is caused by traps. In a finite system there is a given maximal depth of the traps (maximal waiting time) truncating their power-law waiting time distribution in such a way that beyond the maximal waiting time the diffusive behavior of the complex system tends to normal. Note also that the truncation of the waiting-time distribution demonstrates features of weak ergodicity breaking in motion of lipid granules [147]. The above mentioned relaxation functions are only partial cases of the more universal HN law.

6.1 Simple tempered relaxation

The tempered $L\alpha S$ process [148, 149] is characterized by the following Laplace image of its pdf

[TABLE]

where the stability parameter $0<\alpha\leq 1$ and the tempering parameter $\delta\geq 0$ are constants. If $\delta$ equals to zero, the tempered $L\alpha S$ process becomes simply $L\alpha S$. In other words, the parameter $\delta$ provides just a truncation of the ordinary, long-tailed totally skewed $L\alpha S$ distribution. The truncation leads to the random process having all moments finite. Formula (140) describes probabilistic properties of the tempered process in terms of the internal (operational) time. Its inverse process may be used as a subordinator. The pdf $g_{\alpha}(\tau,t)$ of the subordinator depends on the real physical time and describes the first passage over the temporal limit $t$. Its Laplace transform reads

[TABLE]

The inverse tempered $L\alpha S$ process accounts for motion alternating with stops so that the temporal intervals between them are random and with heavy tails in density. The main feature of the process is that it occurs only for short times [142].

Let the parent process $X(\tau)$ have the pdf $h(x,\tau)$. Then the pdf of the subordinated process $X[S(t)]$ obeys the integral relationship between the pdfs of the parent and directing processes, $X(\tau)$ and $S(t)$, respectively,

[TABLE]

The pdf $p(x,t)$ has the most simple form in the Laplace space

[TABLE]

For $\delta=0$ the above expression becomes equal to $u^{\alpha-1}\tilde{h}(x,u^{\alpha})$ and hence corresponding to the pdf without tempering effects.

Let the ordinary Fokker-Planck equation (FPE) $\partial h(x,\tau)/\partial\tau=\hat{L}(x)\,h(x,\tau)$ describe spatio-temporal evolution of a particle subject to the operation time $\tau$. Acting with the operator $\hat{L}(x)$ on the Laplace image $\tilde{p}(x,u)$ in (143), we find

[TABLE]

where $q(x)$ is an initial condition. Using the formal integral representation of the FPE

[TABLE]

and taking the inverse Laplace transform of Eq.(144), we obtain the explicit form of the memory kernel $M(t)$ [142], namely

[TABLE]

For $t\ll 1$ (or $\delta\to 0$) function (146) takes the power form $t^{\alpha}/\Gamma(\alpha)$. However, for $t\gg 1$ (or $\alpha\to 1$) $M(t)$ becomes constant and, as a result, Eq. (145) transforms into the integral form of the ordinary FPE.

To find the characteristics of the tempered diffusion, we consider the Fourier transform of the diffusion process $X[S(t)]$:

[TABLE]

To expose the properties of the “tempered relaxation” function $\phi(t)$, we use the frequency-domain representation (5). Then, for the relaxation under the inverse tempered $L\alpha S$ process, the shape function (5) takes the form

[TABLE]

where $0\leq\sigma\leq 1$ is a constant.

According to Eq. (148), for $\sigma=0$ the shape function describes the CC law. If $\alpha=1$, the function simplifies to the frequency-domain D formula. In the case of $\sigma=1$ it has the CD form. The relaxation directed by the inverse tempered $L\alpha S$ process takes an intermediate place between the superslow relaxation and the exponential one. Such a type of evolution is observed in relaxation experiments (see, for example, [11]).

It follows from experimental investigations [10, 11] that the complex dielectric susceptibility $\chi(\omega)=\chi^{\prime}(\omega)-i\chi^{\prime\prime}(\omega)$ of most dipolar substances demonstrates a peak in the loss component $\chi^{\prime\prime}(\omega)$ at a characteristic frequency $\omega_{p}$, and it is characterized by high- ($\omega\gg\omega_{p}$) and low-frequency ($\omega\ll\omega_{p}$) dependencies (3) and (4), respectively. The tempered relaxation shows

[TABLE]

where $A=1-\sigma^{\alpha}$, $B=(\sigma^{2}+\omega^{2}/\omega_{p}^{2})^{\alpha/2}$ and $C=\alpha\arctan(\omega/\delta)$. For small $\omega$ it is easy to see that $\lim_{\omega\to 0}\chi_{\rm temp}^{\prime}(\omega)\sim\omega$ and $\lim_{\omega\to 0}\chi_{\rm temp}^{\prime\prime}(\omega)\sim 1$ whereas for large $\omega$ we get $\lim_{\omega\to\infty}\chi_{\rm temp}^{\prime}(\omega)\sim\omega^{-\alpha}$ and $\lim_{\omega\to\infty}\chi_{\rm temp}^{\prime\prime}(\omega)\sim\omega^{-\alpha}$. This implies that

[TABLE]

where $n=1-\alpha$, that is in agreement with “energy criterion” hypothesis formulated on the basis of experimental results [10, 11]. However, for small $\omega$ we come to

[TABLE]

This means that the energy lost per cycle does not have a constant relationship to the extra energy that can be stored by a static field. Such an asymptotic behavior suggests (Figure 13) that the tempered relaxation takes an intermediate place between the D, CC and CD types of relaxation.

6.2 Tempered relaxation with clustering patterns

In real physical systems, under influence of the surroundings, the couplings between relaxing entities may differ from “place to place”. To describe this experimental result, we use two effects in the compound subordination of random processes, namely tempering and coupling, suggested in [150]. Tempering of the random subordinator influences on short-time evolution in such a system, and coupling gives an appropriate long-term trend. It should be also noticed that in this scheme the random processes, applied for the subordinator construction, have finite integer moments. Then the relaxation pattern under the tempering and the coupling in the frequency domain reads

[TABLE]

where the stability parameter $0<\alpha\leq 1$, the coupling parameter $0<\gamma\leq 1$ and the tempering parameter $\sigma\geq 0$ are constants. In this case the asymptotic behavior of the susceptibility $\chi(\omega)$ has different (independent) low- and high-frequency power tails:

[TABLE]

In the sense of arbitrary exponents $0<m=\gamma\leq 1$ and $0<(1-n)=\alpha\leq 1$ the analysis of the above-mentioned relaxation pattern shows that the same material (for example, Nylon 610 and Glycerol) can exhibit both the $m<1-n$ and $m>1-n$ relations under different temperature/pressure conditions, and the change between $m>1-n$ and $m<1-n$ is also observed in Polyvinylidene fluoride for different susceptibility peaks (see Table 5.1 in [10]). The form (150) allows one to fit the whole range of the two-power-law spectroscopy data with independent low- and high-frequency fractional exponents being free parameters in the model.

It is also of interest the time-domain description for the relaxation function given by (149). It takes the form

[TABLE]

where $G_{\alpha,\gamma}(\tau,z)$ is the Dirichlet average of the two-parameter Mittag-Leffler function $E_{\alpha,\beta}(x)$ with $\alpha,\beta>0$ [107], namely

[TABLE]

Recall that many special functions of mathematical physics are expressed in terms of an average [151]. On the other hand, it is useful to remark that Eq.(151) is very similar to the relaxation function derived in [119], if one accepts $E_{\alpha,\alpha}[(\delta^{\alpha}-\tau^{-\alpha}_{p})\tau^{\alpha}]$ instead of $G_{\alpha,\gamma}(\tau,z)$. This is clear because the relaxation function of [119] is a particular case of Eq.(151) with $\gamma=1$. In this connection, it should be pointed out that for $\gamma=1$ the Dirichlet average kernel transforms into the Dirac delta-function.

Now we consider what possibilities for fitting of the experimental data gives the shape function (149). Table 3 just serves for this purpose. Firstly, for the corresponding values of the parameters $\alpha$, $\gamma$ and $\sigma$, the relaxation function (149) describes the D, CC and CD types of relaxation. Taking $\gamma=1$, we come to the pseudo CD relaxation derived in [142]. Next, if $\sigma$ equals to zero, then Eq. (149) takes the form (67), termed as JWS in [152]. The special place is accepted by two remained cases called the pseudo CC and the mirror CD. They look like the conventional CC and CD types of relaxation, respectively, but there are also differences which we discuss below.

The most general case of the relaxation fitting function Eq.(149) assumes arbitrary values of its parameters in the framework of their boundaries $0<\alpha\leq 1$, $0<\gamma\leq 1$ and $\sigma\geq 0$. What features does it demonstrate? The loss curves given by this case in a log-log representation will be asymmetric with respect to the loss peak frequency (material constant) $\omega_{p}$, but they will be more flattened than the true CD curve and its mirror image (reflection) with a comparable asymmetry. The same is observed for the CC plot (Figure 14). The function (149) makes it possible to give an expression of $\chi(\omega)$ corresponding with a CC plot [41] making arbitrary angles $\pi\gamma/2$ and $\pi\alpha/2$ with the $\chi^{\prime}(\omega)$ axis at the low-frequency and high-frequency ends, respectively (see more details in [83]). Under $\gamma=1$ and $\sigma=1$ the expression (149) reduces to the CD case, and for $\alpha=1$ we arrive at the mirror CD picture. The classical D relaxation time appears for $\alpha=\gamma=1$. As it is shown in [41], the shape of the CC plot is wholly determined by the intersection angles with the $\chi^{\prime}(\omega)$ axis, and they often are obtained from experimental data of relaxation.

To sum up, the relaxation pattern given by Eq.(149) satisfies the “energy criterion” [10] indicating origins of scaling properties in the evolution of complex stochastic systems, namely $\lim_{\omega\to 0}\chi^{\prime\prime}(\omega)/(\chi^{\prime}(0)-\chi^{\prime}(\omega))=\tan(m\pi/2)$ and $\lim_{\omega\to\infty}\chi^{\prime\prime}(\omega)/\chi^{\prime}(\omega)=\cot(n\pi/2)$. This means that the energy lost per cycle has a constant relationship to the extra energy that can be stored by a static field. However, the D response is behind this property. Interestingly, the CD response and its mirror reflection support the criterion only either on high frequencies or on low ones, respectively, as their behavior on either short or long times is similar to the D response.

7 Subordination impact in understanding of relaxation scenarios

The subordination approach takes a central place in the theory of relaxation based on the stochastic point of view. The main feature of all the dynamical processes in complex systems is their random background. Subordination describes a transformation of one random process to another where the complex system under consideration is just such a converter.

Consider the kinetic equation for a general type of Markovian processes in the form

[TABLE]

where $\hat{\bf W}$ denotes the transition rate operator. Depending on the context [153], the operator $\hat{\bf W}$ can be a differential operator with respect to space coordinates in the case of diffusion or a matrix of transition rates for relaxation phenomena. This equation defines the probability $p(t)$ for the system evolution from one state into others. Next, we determine a new process governed by an inverse infinitely divisible process with the Laplace image (88), namely

[TABLE]

The Laplace transform $\tilde{p}_{\Psi}(s)$ reads

[TABLE]

and leads to

[TABLE]

From this it follows

[TABLE]

where the kernel $M(t)$ is the memory function defined earlier by formula (92). This analysis shows resemblance between the diffusion scenarios and the scheme of two-state relaxation.

To fix the spatial structure of clusters in time, it is useful to pass from the subordinated random processes to the CTRW limit of random variables in one instance. The procedure results in factorization of the space and time characteristics. In what follows, the mixture of symmetric $L\eta S$ (standard normal for $\eta=2$) and completely asymmetric $L\alpha S$ laws has the same distribution as the diffusion front $\tilde{R}(t)=X_{\eta}(S_{\alpha}(t))$, where $X_{\eta}(\tau)$ belongs to the class of symmetric $L\eta S$ processes ($0<\eta\leq 2$), and it is subordinated by the inverse $L\alpha S$ process $S_{\alpha}(t)$ ($0<\alpha\leq 1$). Really, using the 1/$\eta$-self-similarity property of $X_{\eta}(\tau)$, i. e.

[TABLE]

and the $\alpha$-self-similarity of the independent process $S_{\alpha}(t)$ as

[TABLE]

where

[TABLE]

is the $1/\alpha$-self-similarity of the $L\alpha S$ motion $T_{\alpha}(\tau)$, we obtain

[TABLE]

what gives the desired mixture of the symmetric $L\eta S$ and completely asymmetric $L\alpha S$ laws. Here $C_{\alpha},C_{\eta}$ are constants, and $Q_{\alpha,\eta}$ is a random variable connected with the pdf $\vartheta_{\alpha}(x)$ ($0<\alpha\leq 1$), which describes random variables inverse to $L\alpha S$ ones. Observe that the Laplace transform of $\vartheta_{\alpha}(x)$ gives the well-known one-parameter Mittag-Leffler function, i. e.

[TABLE]

According to above relations, the relaxation function takes the CC temporal form

[TABLE]

obtained earlier in [109] for $\eta=2$. The time constant $\tau_{p}=(Ck)^{-\eta/\alpha}$, where $C=(C_{\alpha}C_{\eta})^{1/\eta}$. Note that the random variable $Q_{\alpha,\eta}$ plays the same role (in a generalized sense) as relaxation rates of clusters in the space representation considered above. In this case the cluster structure dominates only on the microscopic level what cannot be said about cooperative regions.

The diffusion front $\tilde{R}(t)=X_{\eta}(Y_{\gamma}[S_{\alpha}(t)])$ is determined by the compound subordination [79, 135], where the process $S_{\alpha}(t)$ subordinates the process $Y_{\gamma}(y)$ independent on $S_{\alpha}(t)$, and the process $Y_{\gamma}(y)=U_{\gamma}[S_{\gamma}(y)]$ itself is expressed in terms of $U_{\gamma}(s)$, being a strictly $L\gamma S$ motion. Hence, we have

[TABLE]

being inverse to $U_{\gamma}(s)$. As the process $U_{\gamma}(y)$ is one-self-similar that we get

[TABLE]

Fixing of the space structure (clusters and super-clusters) in the system at one instance is equivalent to

[TABLE]

The random variable $Y^{O}_{\gamma}(1)\stackrel{{\scriptstyle d}}{{=}}1/{\cal B}_{\gamma}$ is reciprocal to the beta-distributed variable ${\cal B}_{\gamma}$. The beta distribution with parameters $\gamma$ and $1-\gamma$ is known also as the generalized arcsine distribution. In this scenario the HN relaxation function is written as

[TABLE]

where $h_{\gamma}(x)=[\Gamma(\gamma)\Gamma(1-\gamma)]^{-1}x^{-1}(x-1)^{-\gamma}$ is supported for $x>1$ and 0 otherwise. Here the integrand function $E_{\alpha}\left[-(t/\tau_{p})^{\alpha}x\right]$ describes an evolution in the distributions of microscopic quantities (like relaxation rates of clusters) whereas the distribution $h_{\gamma}(x)$ relates to the mesoscopic level of relaxation rates in super-clusters. The analogous analysis may be also provided for the JWS relaxation, but this makes the similar derivations, so we omit them.

In this context it would be useful to recall that the subordinator $U_{\alpha}(t)$ remains the type of relaxation exponential per se [154]. As it has been shown in [155], the stretched exponential relaxation pattern appears from the other subordinator undergoing the one-self-similarity

[TABLE]

In this case the complex system has neither cluster nor super-cluster structures, but all relaxation rates of system entities obey a $L\alpha S$ distribution.

Taking $\gamma=1$, the diffusion front $\widetilde{R}(t)=X_{\eta}(Y_{\gamma}[S_{\alpha}(t)])$ is simplified to $\widetilde{R}(t)=X_{\eta}(S_{\alpha}(t))$. Then Eq.(156) transforms into Eq.(155), and we get a cluster structure in the disordered system without any contribution of super-clusters. In this case the relaxation function is the CC law. On the other hand, for $\alpha=1$ the diffusion front reads $\widetilde{R}(t)=X_{\eta}(Y_{\gamma}(t))$. In this situation the interactions among super-clusters play a dominating role. Their heavy tails of relaxation rates lead to the CD relaxation pattern.

8 Outline

The manifestations of many-body effects in the relaxation of dipolar systems are fundamental and universal, independent of the physical and chemical structures of their interacting entities. The interactions are usually non-trivial and anharmonic, resulting in complex (chaotic) dynamics in the phase space spanned by the coordinates and momenta of the interacting units. The presence of the many-body effects in the dynamics cannot be ignored. This means that the dynamical processes in the such systems have a stochastic background. Nevertheless, the characterization of complex systems wholly is deterministic, in the form of their universal relaxation patterns. The main feature of relaxing complex systems is that their relaxation response is non-exponential in nature. All types of the empirical functions used to fit the relaxation data exhibit the fractional-power dependence of the dielectric responses in frequency and time. It is worth noticing that this unique property is independent on any special details of examined systems. In this framework one can expect that the macroscopic behavior of the complex systems is governed by “averaging principles” like the law of large numbers to be in force. However, the problem is that the macroscopic evolution of complex systems is not attributed to any particular entity taken from those forming the system. The description of the relationship between the local random characteristics of complex systems and the universal deterministic empirical laws of relaxation is performed by the limit theorems of probability theory. Applications of this approach to relaxation and transport for large classes of (physical, chemical, and other) phenomena [10, 30, 11], involving different types of self-similar random processes, has turned out to be very successful in recent years [84, 92, 152, 68, 87, 69, 93]. Especially, very fruitful appeared the idea of subordination description which allows one to characterize a stochastic transport of particular excitation mode in complex systems. The general probabilistic formalisms treat relaxation of the complex systems regardless of the explicit representation of local interactions, the detailed relations between relaxing entities, interaction ranges and the transport efficiencies. In a natural way, they give efficient methods for evaluating the dynamical averages of the relaxation processes. In this article, utilizing the stochastic tools, we have derived the most-known empirical relaxation laws, characterized their parameters, connected the parameters with local random characteristics of the relaxation processes, reconstructed the internal random structure of relaxing systems, justified the energy criterion, demonstrated the transition from analysis of the microscopic random dynamics in the systems to the macroscopic deterministic description by integro-differential equations. It should be noticed that the classical methods of statistical physics take into account the central limit theorem in respect to the probability distributions having finite moments. However, this concept does not help to clarify the nature of relaxation phenomena in complex systems. The above consideration has an evident advantage over the traditional models and goes behind the classical statistical physics. In our approach, to look at the problem of anomalous time evolution of complex systems from different sides, we used several alternative stochastic methods (see Table 4). Such an attempt has some advantages. Firstly, it sets the range of validity of the approaches. Secondly, once this range was established, it allows us to extend the results to (in principle) arbitrary relaxing systems. Thus, joining different methods may lead to a deeper understanding of the complex dynamics found in nature.

Acknowledgements

A.S. is grateful to the Faculty of Pure and Applied Mathematics and the Hugo Steinhaus Center for pleasant hospitality during his visit in Wrocław University of Science and Technology. He is also grateful for a partial support from the NCN Maestro grant No. 2012/06/A/ST1/00258.