Universality and self-similar behaviour of non-equilibrium systems with non-Fickian diffusion

Analytical approaches describing non-Fickian diffusion in complex systems are presented. The corresponding methods are applied to the study of statistical properties of pyramidal islands formation with interacting adsorbate at epitaxial growth. Using the generalized kinetic approach we consider universality, scaling dynamics and fractal properties of pyramidal islands growth. In the framework of generalized kinetics, we propose a theoretical model to examine the numerically obtained data for averaged islands size, the number of islands and the corresponding universal distribution over the island size.


Introduction
It is well known that complex physical, chemical and/or biological systems manifest strong interactions between their elements (atoms, molecules, individuals/organisms) leading to self-organization processes with anomalous diffusion, nonequilibrium phase transitions, patterning and the formation of self-similar objects such as fractals. The cooperative interactions in such systems result in a reduction of a large number of degrees of freedom binding the subunits of these systems by means of self-organization into synergetic entities. These synergetic systems admit low-dimensional description by the corresponding Fokker-Planck equation (FPE). As was shown previously (see for example, [1]), systems manifesting a collective behavior of interacting species are described in terms of nonlinear FPE (NFPE) capturing the essential dynamics underlying the observed phenomena with anomalous transport and generalized (nonextensive) or Q-deformed statistics [2]. The corresponding nonlinearity in NFPE without a drift component, leads to non-Fickian diffusion with density/concetration dependent dispersal resulting in anomalous dynamics as was shown experimentally and theoretically (see, for example, [2][3][4][5][6]).
Nonlinear diffusion equation can be generalized by taking quasi-chemical reactions into account. In such a case, one arrives at nontrivial scenarios for nonequilibrium phase transitions [7], pattern formation [8] and delaying dynamics at phase separation processes [9]. For example, while studying the pattern formation phenomena on surfaces at deposition from a gaseous phase, one describes the corresponding system by reaction-diffusion model with field dependent diffusivity. It was shown that nanosize patterns of adsorbate can emerge due to microscopic interactions of the deposited particles [10][11][12][13][14][15][16]. In [17] it was shown that fractal properties of porous-surface condensates are described by the generalized statistics based on NFPE approach and by the corresponding theory of multifractals. In the process of pyramidal islands growth at molecular beam epitaxy, it was found that a structure of pyramids essentially depends on interactions of the elements forming the pattern described by a concentration dependent diffusion coefficient [18][19][20][21]. While studying the arrangement of point defects in solids at particle irradiation according to the swelling rate theory [22,23], it was found that vacancies can arrange into nanosize clusters due to their interactions described by a nonlinear diffusion flux [24][25][26][27]. This effect can lead to abnormal grain growth dynamics when vacancies segregate on the grain boundaries in a stochastic manner [28]. Nonlinear diffusion was experimentally studied for a large variety of chemical and biological systems [29,30]. In [31][32][33] it was shown that nonlinear diffusion leads to anomalous dynamics. From this non-complete literature overview, it follows that in spite of the diversity of this research area, many phenomena in physical, chemical and biological systems have fundamental physical mechanisms in common. A study of interacting mechanisms of elements of complex systems leading to macroscopic self-organization processes with collective behavior of their species remains an urgent problem in modern statistical physics during the last two decades.
In this paper we initially present a generalized kinetic approach permitting to describe the behavior of physical systems with an interaction of their species. Here, we illustrate a self-similar behavior of the main statistical characteristics for the systems described by nonlinear diffusion equation. Next, studying self-organization processes with the surface pattern formation in a system of interacting adsorbate at molecular beam epitaxy, we use the formalism of nonlinear kinetic approach to describe the scaling properties of pyramidal islands growth and the universality of the system behavior.
The work is organized as follows. In section 2 using the generalized master equation we introduce NFPE and discuss the corresponding anomalous dynamics characterized by nonextensive statistics. Here, we consider the nonlinear diffusion equation and study the properties of a related solution. In section 3 we discuss the physical reasons responsible for the emergence of nonlinear diffusion flux in systems with interacting species. Here, (subsection 3.1) we consider a typical model of reaction-diffusion systems with field dependent diffusivity related to the surface pattern formation at molecular beam epitaxy. The original results illustrating the role of particle interactions and scaling dynamics with universal properties of the system described in the framework of NFPE approach are collected in subsection 3.2. We conclude in section 4.

Generalized kinetic approach and nonlinear diffusion equation
The classical Fickian diffusion assumes that there are no interactions between the random walkers representing individuals in a system. If these walkers interact, the diffusivity or the mobility can change in the presence of walkers and become explicitly dependent functions of the local concentrations of the moving particles. This leads to NFPE as a generalization of ordinary linear FPE where transition probabilities depend on occupation probability densities of initial and arrival states. Following [34], NFPE can be constructed from the generalized master equation written for a probability density function (pdf) p(r, t ) to find a particle in the state r at the time t . This equation has the form where w(r, r ) is the transition rate between two states located at r and at r which depends on the nature of interactions between the particle and the bath; it is a function of starting r and arrival r sites. The factor γ(p, p ) is an arbitrary function of the particle populations of both the initial and the arrival sites. This function satisfies the following conditions: γ(0, p ) = 0 meaning that if the initial site is empty, the transition probability is equal to zero; γ(p, 0) 0 requires that if the arrival site is empty, the transition probability should depend on the population of the initial site. In the case γ(p , p) = p and γ(p, p ) = p, we have the standard liner kinetics when equation (1) reduces to the Chapman-Kolmogorov equation.
In the diffusion limit one can expand slowly varying functions γ in equation (1) and obtain NFPE in the form [34] ∂p ∂t where ∇ ≡ ∂/∂r. The drift and diffusion terms are: f (r) ≡ dr r w(r, r ), 2D(r) ≡ dr r 2 w(r, r ). Here, γ(p) = γ(p, p), and the function κ(p) > 0 is defined through the condition

33004-2
Formally, the functions γ(p) and κ(p) can be represented through densities of the initial a(p) and arrival b(p) states as follows: γ(p) = a(p)b(p), κ(p) = a(p)/b(p) [35]. In stationary equilibrium state, the solution p s = p(r, ∞) of equation (2) takes the form where γ(p) remains an arbitrary function; U 0 takes care of normalization condition p s dr = 1. A nonlinearity of the Fokker-Planck equation is defined through the form of the function κ(p): at ln κ(p) = ln p, we move to the Boltzmann-Gibbs statistics; a Q-deformation of the logarithm ln κ → ln Q κ promotes Tsallis In the simplest case with f (r) = 0 and D = const, one can consider the time dependent solution of the nonlinear diffusion equation as reduced NFPE. Here, by taking γ = κ = p one gets As was shown in [29,36,37], this equation has an exact solution in one dimension at |r | r 0 (t /t 0 ) H : where N 0 is the initial number of particles at the origin. At |r | > r 0 (t /t 0 ) H , one has p(r, t ) = 0. It follows that such a nonlinear diffusion equation admits scaling 〈(r − 〈r 〉) 2 〉 ∝ t 2H , whereas ordinary Brownian diffusion corresponds to H = 1/2 with Q = 1. Scaling regimes in a self-affine phase space characterized by D(r, p) = p 1−Q r ∆ were studied in [5]. An automodel solution of the corresponding nonlinear diffusion equation is characterized by the Hurst exponent H in the form H = 1/(3−Q −∆), whereas the related pdf is well-described by the Tsallis form. Therefore, anomalous dynamics can be controlled either by Tsallis parameter Q or by the exponent ∆.
In the case of the constant drift term [ f (r) = const], the corresponding simulations of the time series, satisfying Tsallis statistics, have shown clusterization of time series manifesting a self-similar regime characterized by fractal properties [6,38]. A study of the stochastic dynamics corresponding to NFPE with a power law dependence γ(p) was done in [39]; here, anomalous dynamics was considered using correlation analysis. The relation between NFPE and Levy-type diffusion was discussed in [36,40]. The application of the generalized statistics followed by NFPE with anomalous diffusion emergent at selforganized criticality regimes was studied in [4], where scaling laws were discussed for avalanche sizes dynamics and the corresponding power-law distributions. Therefore, the nonlinear kinetic approach is applicable to a description of a universality and scaling behavior of complex systems characterized by field dependent diffusivity related to interactions between system elements. In the next section we study the pattern formation processes by considering the reaction diffusion model that describes epitaxial growth with interacting adsorbate. We use the nonlinear kinetic formalism to analytically describe the scaling regimes and universality of probability density functions of pyramidal islands over their sizes.

Universality and scaling regimes of pattern formation at epitaxy
In the previous section we have considered the nonlinear diffusion equation, where density dependent diffusion coefficient was defined through occupation probabilities of the initial a(p) and arrival b(p) states. In this section, while studying the pyramidal islands growth at epitaxy, we illustrate that the corresponding nonlinear diffusivity can be immediately obtained considering the interactions of adsorbed particles. Universality and scaling regimes of the arranged pyramidal islands will be studied with the help of NFPE presented in the previous section.

Phase field model of pyramidal islands growth
Considering a system of interacting adsorbate with free diffusion on a surface and the motion of mobile particles caused by their interactions, one can directly obtain the diffusion flux with a field dependent diffusion coefficient. In such a case, instead of the probability density p, the relevant quantity that can be used is a coverage field ρ (concentration of atoms/moleculas adsorbed by the surface). An evolution of this quantity is described by a reaction-diffusion equation of the is the reaction term including adsorption/desorption and/or additional nonequilibrium processes related to the formation of islands, or oxides; J is the diffusion flux.
At epitaxial growth, the reaction term R(ρ) is defined through a deposition rate characterized by the flux F 0 of arriving particles on a surface and by desorption reaction −ρ/τ, where τ is the relaxation time.
Next, we assume that desorbed particles have lateral interactions described by the attractive potential Assuming that an interaction length r 0 between two interacting species is small compared to a diffusion length = Dτ 0 , i.e., r 0 , one can write U (ρ(r )) − ρ(r ) with = u(r )dr . The presented formalism is general and can be applied to a large class of systems: semiconductors and metals, where the quantity depends on concrete material properties. Using these assumptions one can rewrite the total diffusion flux (7) through the free energy functional F [ρ] in the standard form Therefore, the reaction diffusion model of a system with an interacting adsorbate can be written in the where dimensionless quantities t = t /τ 0 , r = r/ , F = F 0 τ 0 , θ = T /T 0 , ε = /T 0 are used; T 0 is the bath temperature. Next, we drop the primes for convenience. This approach is widely used to study the pattern formation processes in chemical systems, at condensation under deterministic (see [10][11][12]) and stochastic conditions (see [13][14][15][16]), at epitaxial growth [20,21], at the formation of point defect clusters due to irradiation effect [24][25][26]41] and at other physical systems manifesting interactions between their elements. The obtained model cannot be used immediately to model the pyramidal islands formation because it does not take into account discrete steps (sharp interfaces) between terraces in pyramids. The problem of sharp interface modelling at the pyramidal islands formation can be solved using the phase field approach proposed by Liu and Metiu [19] and Karma and Plapp [18]. The idea of the phase field approach lies in introducting an order parameter φ(r, t ) that indicates the phase at a particular position. In our case, the phase field φ describes the surface height in units of monoatomar layers. According to [19], local stable minima of the order parameter relate to terraces whereas rapid spatial variation of φ corresponds to the positions of steps. In the framework of the phase field approach [18], one has τ φ ∂ t φ = − δH δφ .
Here, τ φ is the characteristic scale of the time of attachment of adatoms at the step. H is the effective Hamiltonian, H = dr[ 2 (∇φ) 2 Here, stands for the width of the step, λ is the dimensionless coupling constant, φ s /2 is the height of the initial substrate. This model reduces to the solid-on-solid model when the coupling between two fields is a constant supersaturation. The term incorporating 1 + cos(π[φ − φ s ]) admits that minima of the free energy H are possible at φ − φ s = 2n + 1, independently of the adatom concentration [18]. According to the phase field approach, the corresponding dynamics of the coverage field is [18]. This model is relevant to the case of a constant adsorbate temperature. It is known that the temperature of the growing surface can be changed locally at adsorption/desorption processes: when atom becomes adatom, the temperature increases locally, it decreases when the desorption of adatom occurs. Moreover, the temperature can be increased due to the effect of the source of atoms described by F . Using the above mechanisms for the temperature variations, one can write the evolution equation for the temperature field in the dimensionless form [21] µ∂ t θ = 1 − θ + χ∆θ + r a F ρ + r b ∂ t ρ. (9) The relaxation of the temperature to the bath temperature is described by two first terms in the righthand side with a relaxation time τ T , where µ = τ T /τ 0 ; χ plays the role of thermal diffusivity. The fourth term describes re-heat of the surface with an intensity r a due to energy exchange with the environment by deposition flux F . This is a standard assumption and it is widely used considering temperature instabilities in chemical reactions (see, for example [41,42]). Such temperature instability can be caused by reorganization of the surface due to annihilation of defects on the surface and their motion to sinks. The fourth term in equation (9) does not directly take into account a change of the temperature for curved steps. Such temperature change can be described by the phase field. In our consideration, this effect is effectively taken into account through the introduction of the last term in equation (9). It also relates to the local heating (∂ t ρ > 0) or cooling (∂ t ρ < 0) during adsorption/desorption processes with intensity r b ; generally, it is responsible for the formation of a curved step.
The total system of three equations describing epitaxial growth generalized by introducing fluctuation terms is [21] where ϑ = τ φ /τ 0 . The last terms in equation (10) are stochastic sources responsible for statistical description of the system dynamics: Using the linear stability analysis, one can find critical values of main control parameters F and ε bounding a domain of their values corresponding to the formation of pyramidal islands (see figure 1). At ε < ε c and F < F c , the averaged phase field 〈φ〉 does not grow with time meaning the formation of a surface with Gaussian profile related to the effect of fluctuation terms in equation (10). In the opposite case, one has ∂ t 〈φ〉 > 0, resulting in an increase of the surface height and in the formation of well organized pyramidal structures on the growing surface. In [20,21], it was shown that the morphology of pyramids essentially depends on the effect of the nonlinear diffusion term in the flux of the adsorbate. Snapshots of

33004-5
a typical evolution of the system are shown in figure 2 at F > F c . In our simulations, we use a square lattice with linear size L = 256 sites and periodic boundary conditions. We choose a time step ∆t = 2.5×10 −4 and the mesh size l = 1. The Gaussian distributed field φ(r, 0) with 〈φ(r, 0)〉 = 0 and 〈[δφ(r, 0)] 2 〉 = 0.1 was taken as initial conditions. For initial coverage and temperature fields we use ρ(r, 0) = 0 and θ(r, 0) = 1.
It is seen that at the initial stages, small islands of adsorbate emerge. Here, the temperature field can be locally changed. During the system evolution, such islands become centers of pyramidal patterns, and the corresponding pyramids connect each other by terraces of equivalent heights (see figure 2).

Results and discussions
Let us study the universal dynamics of the islands growth. To this end, we compute the average islands area 〈s(t )〉 at a half-height of the whole system of pyramids emergent at different times and compute the number of the corresponding islands N (t ). It was found that there is a strong relation between 〈s(t )〉 and N (t ):

33004-6
pdf of the island area distribution N (s, t ) should satisfy the following two criteria: and ∞ 0 sN (s, t )ds = S 0 ≡ const. The first one defines the number of islands, whereas the second one corresponds to the surface conservation law. This law is applicable only at the first stages of the pyramids growth, whereas at late stages one gets only one large constantly growing pyramid. The corresponding dependencies of 〈s(t )〉 and N (t ) are shown in figures 3 (b), (c). It is seen that the scaling exponent β depends on the main system parameters reduced to F and ε, and takes the values smaller or larger than 1 related to the normal growth regime: 〈s(t )〉 ∝ t .  To characterize the distribution of the islands area, we compute the integral distribution function P (s, t ) at first and get the corresponding pdf N = dP /ds using numerical differentiation (see figure 4). From the obtained dependencies it follows that the universal behavior of both P (s, t ) and N (s, t ) at different times is realized. Moreover, the corresponding fitting procedure gives logarithmic approximation of P (s) and algebraic dependencies in the form of the Zipf law N (s) = C 2 s −1 at different times, where the normalization parameter is a function of time, i.e., C 2 = C 2 (t ).
To describe the obtained universal dynamics and the universal behavior of the distribution function over islands area that satisfy the above area conservation condition, we apply the formalism of NFPE to our system. Since the island area grows diffusively (by attachment/deattachment interacting atoms), we start with the nonlinear diffusion equation for N (s, t ) written in the form where we assume the form D(N ) = D 0 N 1−Q ; D 0 = const, ∇ s ≡ ∂/∂s for generalized diffusivity. To obtain the corresponding time dependent solution, we consider the system in an automodel regime assuming: is the scaling function that measures the diffusion package smearing. Using the normalization condition and the area conservation law, we immediately get = −2.
Substituting such constructions into equation (11), we can separate a time dependent part and a part depending on y: Here, λ 0 is the constant related to the time dependence. From the first equation we get a relation between Tsallis parameter Q and the scaling exponent β in the form To In further consideration, the quantity 1/D 0 λ 0 can be considered as a small parameter of the theory. It allows us to obtain a reduced pdf of the form ϕ ∝ y −1 with C 0 = 0 and δ = Q − 1. In such a case, the desired pdf is where t 0 and s 0 depend on D 0 and relate to a normalization condition. It is seen that depending on the parameter Q, one gets different dynamics for islands area growth whereas the corresponding area distribution remains universal, independently of the Tsallis parameter Q. According to the obtained numerical data [see figures 3 (b), (c)] and equation (13) for the Tsallis parameter, one has 1 Q < 2. For the normal islands growth characterized by β = 1, we have Q = 3/2. Therefore, the delayed dynamics observed at an elevated deposition rate F and at interaction energy ε is characterized by Q ∈ [1, 3/2), whereas the enhanced dynamics observed at low F or ε relates to Q ∈ (3/2, 2). Using the above formalism for nonlinear diffusion equation for N (s, t ) we can write down the Langevin equation corresponding to equation (11) following [5,39]. It takes the form where ξ(t ) is the white noise having standard properties. Using its formal solution together with the correlator 〈ξ(t )ξ(t )〉 = δ(t − t ) in the automodel regime, we obtain 〈[s(t ) − s(0)] 2 〉 ∝ t 2β(Q−1)+1 . Comparing this time dependence with a priori known 〈[s(t ) − s(0)] 2 〉 ∝ t 2β , one immediately arrives at the relation (13). Testing numerical solutions of the Langevin equation (16) we have found a good agreement between the analytical results and numerical data for time dependencies of 〈s 2 (t )〉 ∝ t β n , where β n is the fitting exponent, and the corresponding pdfs with hyperbolic approximation s −1 at fixed times (see figure 5).
In our simulations, we took fixed Q and obtained β(Q) β n (Q) with errors up to 5% [see protocols in figure 5 (a)]. The corresponding data for pdfs shown in figure 5 (b) allow one to elucidate that there is a domain for s where universal behavior of pdf is realized at different times and at different exponent Q. It follows that with an increase in Q one gets a growing interval for s/s 0 where universal Q-independent behavior of pdf is realized.
The jagging of time series s(t ) can be studied using multifractal detrended fluctuation analysis (see [43]) as a generalization of the standard multifractal theory [38]. Following [43] we can obtain a set of fractal exponents h(q) and a singularity spectrum f (α) = q[α − h(q)] + 1 for multifractals, where α = h(q) + qh (q) is the singularity strength, q is the index for multifractality. It is known that the quantity h(q = 2) coincides with the Hurst exponent 0 H 1 measuring the smoothness of the time series [44] and defining the fractal dimension of time series as D f = 2 − H [38].

Conclusions
Using a nonlinear kinetic approach we have described universality and scaling properties of pyramidal islands formation at epitaxial growth. In the framework of the phase-field model for pyramidal islands growth in systems with interacting adsorbate, we have analyzed the dynamics of the number of islands, their size and the corresponding distribution functions. We proposed a generalized theoretical model for the island size dynamics manifesting a self-similar behavior and universality of the probability density over the island size. The corresponding dynamics is studied using multifractal time series analysis. It is shown that multifractality of the corresponding time series is related to time correlations.