Influence of external potentials on heterogeneous diffusion processes
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug Influence of external potentials on heterogeneous diffusion processes
Rytis Kazakeviˇcius ∗ and Julius Ruseckas Institute of Theoretical Physics and Astronomy, Vilnius University, Saul˙etekio 3, LT-10222 Vilnius, Lithuania
In this paper we consider heterogeneous diffusion processes with the power-law dependence ofthe diffusion coefficient on the position and investigate the influence of external forces on the re-sulting anomalous diffusion. The heterogeneous diffusion processes can yield subdiffusion as wellas superdiffusion, depending on the behavior of the diffusion coefficient. We assume that not onlythe diffusion coefficient but also the external force has a power-law dependence on the position. Weobtain analytic expressions for the transition probability in two cases: when the power-law exponentin the external force is equal to 2 η −
1, where 2 η is the power-law exponent in the dependence ofthe diffusion coefficient on the position, and when the external force has a linear dependence onthe position. We found that the power-law exponent in the dependence of the mean square dis-placement on time does not depend on the external force, this force changes only the anomalousdiffusion coefficient. In addition, the external force having the power-law exponent different from2 η − PACS numbers: 05.40.-a, 02.50.-r, 05.10.Gg
I. INTRODUCTION
There are many systems and processes where the time dependence of the centered second moment is not linearas in the classical Brownian motion. Such family of processes is called anomalous diffusion. In one dimension theanomalous diffusion is characterized by the power-law time dependence of the mean square displacement (MSD) [1] h (∆ x ) i = K α t α . (1)Here K α is the anomalous diffusion coefficient. When α = 1, this time dependence deviates from the linear functionof time characteristic for the Brownian motion. If α <
1, the phenomenon is called subdiffusion. Occurrence ofsubdiffusion has been experimentally observed, for example, in the behavior of individual colloidal particles in two-dimensional random potential energy landscapes [2]. It has been theoretically shown that active particles moving atconstant speed in a heterogeneous two-dimensional space experience self trapping leading to subdiffusion [3]. Usuallyit is assumed that subdiffusive behavior is caused by the particle being trapped in local minima for prolonged timesbefore it escapes to a neighboring minima. Continuous time random walks (CTRWs) with on-site waiting-timedistributions falling slowly as t − α − predict a subdiffusive behavior [4, 5].Superdiffusion processes, characterized by the nonlinear dependence (1) of the MSD on time with the power-lawexponent in the range 1 < α <
2, constitute another subclass of anomalous diffusion processes. Superdiffusion isobserved, for example, in vibrated dense granular media [6]. Theoretical models suggest that supperdiffusion canbe caused by L´evy flights [4]. L´evy flights resulting in a superdiffusion can be modeled by fractional Fokker-Planckequations [7] (or Langevin equations with an additive L´evy stable noise). In many experimental studies it is onlypossible to show that signal intensity distribution has L´evy law-tails: distribution function of turbulent magnetizedplasma emitters [8] and step-size distribution of photons in hot vapors of atoms [9] have L´evy tails. This indirectlyshows that in these systems superdiffusion could be found.Anomalous diffusion does not uniquely indicate the processes occurring in the system, because there are differentstochastic processes sharing the behavior of the MSD (1). The physical mechanisms leading to the deviations fromthe linear time dependence of the MSD can depend on the system or on the temporal and spatial ranges underconsideration. For example, diffusion described by CTRW has been observed for sub-micron tracers in biologicalcells [10–12], structured coloidal systems [13] and for charge carrier motion in amorphous semiconductors [5, 14].Fractional Brownian motion and fractional Langevin equations has been used to model the dynamics in membranes[15, 16], motion of polymers in cells [17], tracer motion in complex liquids [18, 19]. Diffusion of even smaller tracersin biological cells has been described by a spatially varying diffusion coefficient [20]. ∗ [email protected] Recently, in Refs. [21–24] it was suggested that the anomalous diffusion can be a result of heterogeneous diffusionprocess (HDP), where the diffusion coefficient depends on the position. Spatially dependent diffusion can occur inheterogeneous systems. For example, heterogeneous medium with steep gradients of the diffusivity can be created inthermophoresis experiments using a local variation of the temperature [25, 26]. Mesoscopic description of transportin heterogeneous porous media in terms of space dependent diffusion coefficients is used in hydrology [27, 28]. Inturbulent media the Richardson diffusion has been described by heterogeneous diffusion processes [29]. Power-lawdependence of the diffusion coefficient on the position has been proposed to model diffusion of a particle on randomfractals [30, 31]. In bacterial and eukaryotic cells the local cytoplasmic diffusivity has been demonstrated to beheterogeneous [20, 32]. Motion of a Brownian particle in an environment with a position dependent temperature hasbeen investigated in Ref. [33]. In random walk description the spatially varying diffusion coefficient can be includedvia position dependence of the waiting time for a jump event [34], the position dependence occurs because in theheterogeneous medium the properties of a trap can reflect the medium structure. This is the case for diffusion onfractals and multifractals [35]. Inhomogeneous versions of continuous time random walk models for water permeationin porous ground layers were proposed in Ref. [36]. Heterogeneous diffusion process might be applicable to describeanomalous diffusion in such systems.The goal of this paper is to consider HDPs with the power-law dependence of the diffusion coefficient on theposition and to analytically investigate the influence of external forces on the resulting anomalous diffusion. Theinfluence of the external forces on HDPs has not been systematically analyzed. In this paper we assume that theforces are characterized by a power-law dependence on the position. Such forces can arise in various systems: Inmany cases the potentials causing the deterministic forces are power-law functions of position, for example linearand harmonic potentials. Power-law potential with arbitrary power-law exponent acting on a nanoparticle have beencreated in Ref. [37]. Logarithmic potentials yielding forces behaving as x − have been applied to describe dynamicsof particles near a long, charged polymer [38], momentum diffusion in dissipative optical lattices [39], long-rangeinteracting systems [40], bubbles in DNA molecules [41]. As we demonstrate, the external force having a certain valueof the power-law exponent does not restrict the region of diffusion and change only the anomalous diffusion coefficientwithout changing the scaling exponent α . Other values of the power-law exponent in the external force can lead tothe exponential cut-off of the probability density function (PDF) and restrict the region of diffusion, limiting the timeinterval when the anomalous diffusion occurs. We expect that the results obtained in this paper may be relevant fora more complete understanding of anomalous diffusion processes.The paper is organized as follows: In Sec. II we summarize the main properties of HDPs with the power-lawdependence of the diffusion coefficient on the position. The influence of an external force with a particular value ofthe power-law exponent is investigated in Sec. III. This value of the power-law exponent is the same as in the driftcorrection for transformation from the Stratonovich to the It´o stochastic equation. We consider external forces havingother values of the power-law exponent in Sec. IV. Sec. V summarizes our findings. II. FREE HETEROGENEOUS DIFFUSION PROCESS
Heterogeneous diffusion process with the power-law dependence of the diffusion coefficient on the position, intro-duced in Ref. [21], is described by the Langevin equation dx = σ | x | η ◦ dW t . (2)Here η is the power-law exponent of multiplicative noise, σ is the amplitude of noise and W t is a standard Wienerprocess (Brownian motion). This stochastic differential equation (SDE) is interpreted in Stratonovich sense. Formathematical convenience and for further generalization we transform Eq. (2) to the Itˆo convention: dx = 12 σ η | x | η − xdt + σ | x | η dW t . (3)The Fokker-Planck equation corresponding to the SDE (2) is [42] ∂∂t P ( x, t ) = 12 σ ∂∂x (cid:20) | x | η ∂∂x ( | x | η P ( x, t )) (cid:21) . (4)With the reflective boundaries at small positive x = x min and large x = x max , Eq. (4) leads to the steady-state PDF P ( x ) ∼ x − η . Without such boundaries the time-dependent solution of Eq. (4) with the initial condition P ( x,
0) = δ ( x )is given by a stretched (when η >
0) or compressed (when η <
0) Gaussian [21] P ( x, t ) = | x | − η √ πσ t exp (cid:18) − | x | − η ) − η ) σ t (cid:19) . (5)When η <
1, this solution describes exponential cut-off at large values of x and power-law behavior at small valuesof x . The position of the cut-off moves towards large values of x with increase of time t . In contrast, when η > x and power-law behavior at large values of x . Theposition of cut-off moves towards smaller values of x with increase of time t .In Ref. [21] it has been demonstrated that Eq. (2) (or, equivalently, Eq. (3)) leads to the power-law time dependenceof the MSD h x ( t ) i ∼ ( σ t ) − η . (6)This behavior of the MSD means that HDP described by Eq. (2) yields superdiffusion for 1 > η > η <
0. For η > x , leading to a progressive acceleration of the diffusingparticle. In contrast, the diffusivity decreases with increasing x when η <
0. For η > η = 1, the MSD grows not as a power-law of time, but exponentially [43, 44].The HDP (2) displays weak non-ergodicity, that is the scaling of time and ensemble averages is different. Specifically,in Ref. [21] it has been shown that the average over the trajectories D δ (∆) E = 1 N N X i =1 δ i (∆) (7)of the the time-averaged MSD δ (∆) = 1 T − ∆ Z T − ∆0 [ x ( t + ∆) − x ( t )] dt (8)scales as D δ (∆) E ∼ ∆ T ηη − . (9)Thus time-averaged MSD depends on the time difference ∆ linearly, in contrast to the power-law behavior of MSDin Eq. (6).Another interesting property of HDPs is the behavior of the distribution of the time-averaged MSD δ of individualrealizations. When η <
0, the distribution of δ decays to zero at δ = 0 [24]. This behavior of the distribution inHDPs is different than the behavior in CTRWs, where there is a finite fraction of immobile particles resulting in thefinite value of the distribution at δ = 0. This difference allows to distinguish between different origins of anomalousdiffusion. III. EXTERNAL FORCE THAT DOES NOT LIMIT THE ANOMALOUS DIFFUSION
In general, not only a random force leading to the diffusion but also a deterministic drift force can be present in theLangevin equation. Therefore, the question arises how external force influences the anomalous diffusion described byheterogeneous diffusion process. To investigate this question we will consider Eq. (3) with an additional drift term.As the results obtained in this Section show, for the appearance of the anomalous diffusion it is sufficient to consideronly positive values of x . Therefore from now on we will write just x instead of the absolute value | x | , assuming thepresence of the boundary that does not allow for the diffusing particle to enter the region x x as the drift correction for the transformationfrom the Stratonovich to the It´o stochastic equation. Such a drift term can arise not only due to an external forcebut can also represent a noise-induced drift [45]. Thus we will generalize the SDE (3) to take the form dx = σ (cid:16) η − ν (cid:17) x η − dt + σx η dW t . (10)Here ν is a new parameter describing the additional drift term. The meaning of the parameter ν is as follows: when thereflective boundaries at small positive x = x min and large x = x max are present, the steady-state PDF is a power-lawfunction of position with the power-law exponent ν , P ( x ) ∼ x − ν . Comparison of Eq. (10) with Eq. (3) shows thatthe free HDP is obtained when ν = η .The SDE (10) has been proposed in Refs. [46, 47] to generate signals having 1 /f β noise in a wide range of frequencies.Such nonlinear SDEs have been applied to describe signals in socio-economical systems [48, 49] and as a model ofneuronal firing [50]. According to Ref. [47], the power-law exponent β in the power spectral density S ( f ) ∼ f − β isrelated to the parameters of SDE (10) as β = 1 + ν − η − . (11)For some values of the parameter ν the SDE (10) can be obtained from Eq. (2) by changing the prescription forcalculating stochastic integrals. Let us consider the SDE dx = σx η ◦ γ dW t , (12)together with the interpretation of stochastic integrals as [45, 51] Z T f ( x ( t )) ◦ γ dW t = lim N →∞ N − X n =0 f ( x ( t n ))∆ W t n , t n = n + γN T . (13)The parameter γ , 0 γ
1, defines the prescription for calculating stochastic integrals. Commonly used values of theparameter γ are γ = 0 corresponding to pre-point Itˆo convention, γ = 1 / γ = 1 corresponding to post-point H¨anggi-Klimontovich [52, 53], kinetic or isothermal convention[54–56]. The integration convention should be determined from the experimental data or derived from another model[57]. The SDE (12) also has been investigated in Ref. [58]. Transformation of Eq. (12) to the Itˆo equation yields[51, 59, 60] dx = σ γ ( η − x η − dt + σx η dW t . (14)This equation has the form of SDE (10) with the parameter ν being ν = 2[ γ + (1 − γ ) η ] . (15)Since 0 γ
1, the range of possible values of the parameter ν obtained by changing the prescription in Eq. (2) islimited: 2 η ν η < ν η when η >
1. In this Section we do not place such restriction on thepossible values of ν . The values of ν outside of this range can be obtained due to the action of an external force.The Fokker-Planck equation corresponding to the SDE (10) is [42] ∂∂t P ( x, t ) = σ (cid:16) ν − η (cid:17) ∂∂x [ x η − P ( x, t )] + σ ∂ ∂x [ x η P ( x, t )] . (16)The time-dependent PDF of the process given by Eq. (10) can be obtained as follows: Transformation of the variable x to a new variable y = x − η (assuming that η = 1) leads to the SDE dy = − σ ′ ν ′ y dt + σ ′ dW t , (17)where ν ′ = η − νη − , σ ′ = | η − | σ . (18)Equation (17) has the form of a Bessel process [61]. This connection with the Bessel process has also been pointedout in Ref. [58]. The known analytic form of the solution of the Fokker-Planck equation ∂∂t P y = 12 σ ′ ν ′ ∂∂y y − P y + 12 σ ′ ∂ ∂y P y (19)corresponding to SDE (17) is [61–63] P ( y, t | y ,
0) = y − ν ′ y ν ′ σ ′ t exp (cid:18) − y + y σ ′ t (cid:19) I − ν ′ +12 (cid:16) yy σ ′ t (cid:17) . (20)Here I n ( z ) is the modified Bessel function of the first kind. This PDF obeys the initial condition P ( y, t = 0 | y ,
0) = δ ( y − y ). For completeness, one possible way of obtaining this solution of Eq. (19) is presented in Appendix A.The PDF (20) can be normalized and represents an acceptable solution only when − ν ′ +12 > −
1, that is when ν ′ <
1. When this inequality is not satisfied, the Bessel process leads to a total absorption at the origin in afinite time [64]. Transforming the variables back we obtain the time-dependent PDF satisfying the initial condition P ( x, t = 0 | x ,
0) = δ ( x − x ): P ( x, t | x ,
0) = x − η − ν x − η + ν | η − | σ t exp − x − η ) + x − η )0 η − σ t ! I ν +1 − η η − x (1 − η ) x (1 − η )0 ( η − σ t ! . (21)Similarly as Eq. (20), the solution (21) is valid when ν +1 − η η − > −
1. This condition is equivalent to ν > η > ν < η <
1. Similar expression of the PDF (21) has been obtained in Ref. [58] by considering the HDP withvarious values of the prescription parameter γ .Using Eq. (21) we can calculate the time-dependent average of a power of x : h x a i x = Z ∞ x a P ( x, t | x , dx = Γ (cid:16) ν − − a η − (cid:17) Γ (cid:16) ν − η − (cid:17) (cid:0) η − σ t (cid:1) a − η ) F a η −
1) ; ν − η −
1) ; − x − η )0 η − σ t ! . (22)Here F ( a ; b ; z ) is the Kummer confluent hypergeometric function. The integral is finite when: i) η > ν > a with a > ν > a <
0; ii) η < ν < a with a < ν < a >
0. Eq. (22) for the a -th momentof x has been derived also in Ref. [58]. In particular, the average of x is equal to h x i x = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) (cid:0) η − σ t (cid:1) − η ) F η −
1) ; ν − η −
1) ; − x − η )0 η − σ t ! (23)and is finite when ν > η > ν < η <
1. The average of the square of x is equal to h x i x = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) (cid:0) η − σ t (cid:1) − η F η −
1) ; ν − η −
1) ; − x − η )0 η − σ t ! (24)and is finite when ν > η > ν < η <
1. When time t is large, that is when x − η )0 η − σ t ≪ , (25)the hypergeometric function in Eq. (22) is approximately equal to 1. Thus for large time the average h x a i x does notdepend on the initial position x and is a power-law function of time: h x a i x ≈ Γ (cid:16) ν − − a η − (cid:17) Γ (cid:16) ν − η − (cid:17) (cid:0) η − σ t (cid:1) a − η ) . (26)As a consequence we obtain that for large time t satisfying the condition (25) the average of the square of x dependson time as t / (1 − η ) . In addition, using Eqs. (23) and (24) we get that the the variance h ( x −h x i ) i = h x i−h x i has thesame dependence on time: h ( x − h x i ) i ∼ t / (1 − η ) . We can conclude, that not only the original HDP equation (2), butalso the equation (10) with additional power-law force exhibit anomalous diffusion. The power-law exponent in thetime dependence is the same as in Eq. (6), the external force considered in this Section does not change the characterof the anomalous diffusion as long as the condition (25) is satisfied. The power-law exponent 1 / (1 − η ) depends only onthe random force in the SDE (10) and does not depend on the parameter ν characterizing the deterministic externalforce. The external force changes only the anomalous diffusion coefficient.We check the analytic results obtained in this Section by comparing them with numerical simulations. The power-law form of the coefficients in SDE (10) allows us to introduce an operational time τ in addition to the physical time t -1 -1 ⟨ x ⟩ t (a) -2 -1 -1 ⟨ x ⟩ t (b) -3 -2 -1 -1 ⟨ x ⟩ t (c) -2 -1 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (d) -2 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (e) -4 -2 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (f) FIG. 1. Dependence of the mean (a,b,c) and variance (d,e,f) on time for various values of the parameters η and ν when theposition of the diffusing particle changes according to Eq. (10). Solid gray lines show numerical result, dashed black lines arecalculated using Eqs. (23) and (24), dotted (red) lines show the power-law dependence on time ∼ t / [2(1 − η )] for (a,b,c) and ∼ t / (1 − η ) for (d,e,f). The parameters are σ = 1 and η = − , ν = − η = , ν = 0 for (b,c); η = , ν = 5 for (c,f).The initial position is x = 1. so that the diffusion coefficient in the operational time becomes constant [65]. The relation between the physical time t and the operational time τ is specified by the equation dt = σ − x − η dτ . For the numerical solution of SDE (10) weuse the Euler-Maruyama scheme with a variable time step ∆ t k = ∆ τ / ( σ x ηk ) which is equivalent to the introductionof the operational time [65]. Thus the numerical method of solution of SDE (10) is given by the equations x k +1 = x k + (cid:16) η − ν (cid:17) ∆ τx k + √ ∆ τ ε k , (27) t k +1 = t k + ∆ τσ x ηk . (28)Here ∆ τ ≪ ε k are normally distributed uncorrelated random variableswith a zero expectation and unit variance. To avoid the divergence of the diffusion and drift coefficients at x = 0in the numerical simulation, we insert a reflective boundary at x = 10 − . This modification is analogous to theregularization of the diffusivity at x = 0, performed in Refs. [21, 24]. When η > ν > x increasesfor small values of x . In this case the operational time introduced in such a way that the coefficient before noisebecomes proportional to the first power of x can lead to a more efficient numerical method. Such an operational timeis introduced by the equation dt = σ − x − η − dτ and the numerical method of solution becomes x k +1 = x k h τ (cid:16) η − ν (cid:17) + √ ∆ τ ε k i , (29) t k +1 = t k + ∆ τσ x η − k . (30)We have calculated the mean and the variance by averaging over 10 trajectories. Comparison of the analyticexpressions (23), (24) with numerically obtained time-dependent mean and variance is shown in Fig. 1. For numericalsimulation we have chosen three different values of the exponent η : η = − , η = , and η = corresponding,respectively, to subdiffusion, superdiffusion and to the localization of the particle. For each case we have chosen avalue of the parameter ν that differs from η of the free HDP. For all numerical simulations the initial position is x = 1. We see a good agreement of the numerical results with analytic expressions. As we can see in Fig. 1, the timedependence of the mean and the variance becomes a power-law for large times and the initial position is forgotten (for -5 -4 -3 -2 -1 -2 -1 P ( x ) x (a) -5 -4 -3 -2 -1 -2 -1 P ( x ) x (b) -4 -2 -2 -1 P ( x ) x (c) FIG. 2. Time-dependent PDF P ( x, t | x ,
0) corresponding to times t = 1 (light gray) and t = 10 (dark gray) for various valuesof the parameters η and ν when the position of the diffusing particle changes according to Eq. (10). Dashed black lines arecalculated using Eq. (21). The dotted line shows the slope x − ν . The parameters are σ = 1 and (a) η = − , ν = −
1; (b) η = , ν = 0; (c) η = , ν = 5. The initial position is x = 1. the parameters used in Fig. 1(b) and (e) the difference between the exact solution and the power-law approximationremains constant, but the relative difference is decreasing).Comparison of the analytic expression (21) for the time-dependent PDF P ( x, t | x ,
0) with the results of numericalsimulation is shown in Fig. 2. To illustrate how the PDF changes with time, the PDF is shown at two different timemoments t = 1 and t = 10. We see a good agreement of the numerical results with the analytic expression for bothtime moments. With increasing time the PDF shifts to the larger values of x when η < x when η > IV. EXTERNAL FORCE LEADING TO AN EXPONENTIAL RESTRICTION OF THE DIFFUSION
Now let us consider the external deterministic force having a power-law dependence on x but with the power-lawexponent different than 2 η −
1. When such a force is positive if the power-law exponent is smaller than 2 η − η −
1, the SDE describing the HDP can be written as dx = σ (cid:18)(cid:16) η − ν (cid:17) x η − + m x m min x η − − m − m x m max x η − m (cid:19) dt + σx η dW t . (31)Here m , m > x min , x max are the parameters of the external force. In Eq. (31) we included three terms in thedrift. Each term has power-law dependence on position x , but with different power-law exponents: equal to 2 η − η − η −
1. The steady-state PDF corresponding to Eq. (31) is P ( x ) ∼ x − ν exp (cid:26) − (cid:16) x min x (cid:17) m − (cid:18) xx max (cid:19) m (cid:27) . (32)We see that the additional terms in Eq. (31) lead to exponential restriction of the diffusion. When x min ≪ x ≪ x max ,the steady-state PDF has the power-law form P ( x ) ∼ x − ν . Confined HDP has been investigated in Refs. [66, 67].Analysis of confined HDP can be relevant for the description of the tracer particles moving in the confinement ofcellular compartments or for the particle traced with optical tweezers that exert a restoring force on the particle [12].We can mathematically obtain one of the terms in Eq. (31) leading to the exponential restriction of the diffusionby transforming the time in the initial equation. Indeed, if we start with the Itˆo SDE dx = a ( x ) dt + b ( x ) dW t (33)and introduce a new stochastic variable z ( t ) = g ( t ) x ( c ( t )) then the SDE for the stochastic variable z becomes dz = (cid:18) dg ( t ) dt zg ( t ) + g ( t ) dc ( t ) dt a (cid:18) zg ( t ) (cid:19)(cid:19) dt + g ( t ) r dc ( t ) dt b (cid:18) zg ( t ) (cid:19) dW t . (34)In Eq. (34) a new term that is proportional to z and to the derivative of g ( t ) appears in the drift. The PDF of thestochastic variable z is related to the PDF of x via the equation P z ( z, t ) = 1 g ( t ) P x (cid:18) zg ( t ) , c ( t ) (cid:19) . (35)Thus, to introduce a new term into Eq. (10), let us start with a stochastic variable y ( t ) obeying the SDE (10) andconsider a new stochastic variable x ( t ) = e µt y (cid:18) κ ( e κt − (cid:19) . (36)Here the functions g ( t ) and c ( t ) are g ( t ) = e µt and c ( t ) = κ − ( e κt − x is dx = (cid:16) µx + σ (cid:16) η − ν (cid:17) e κt − µ ( η − t x η − (cid:17) dt + σe κt − µ ( η − t x η dW t . (37)We can obtain an equation with time-independent coefficients by requiring that κ = 2 µ ( η − . (38)Using this value for the parameter κ the SDE for x becomes dx = (cid:16) µx + σ (cid:16) η − ν (cid:17) x η − (cid:17) dt + σx η dW t . (39)When µ has the same sign as η −
1, this SDE can be written in the form similar to Eq. (31): dx = σ (cid:18) η − ν η − (cid:16) x m x (cid:17) η − (cid:19) x η − dt + σx η dW , (40)where the parameter x m is defined by the equation µ = σ ( η − x η − . (41)Comparing Eq. (40) and Eq. (31) we see that the time transformation considered in this Section introduces anexponential restriction of the diffusion at small values of x when η > x when η < x obeying SDE (40): P ( x, t | x ,
0) = 2 | η − | x η − − e − µ ( η − t x − ν − η x ν − η e ν − η µt × exp − x η − − e − η − µt (cid:16) x − η ) + x − η )0 e − η − µt (cid:17)! × I ν − η η − x η − x (1 − η ) x (1 − η )0 sinh (( η − µt ) ! (42)The conditions of validity of this expression is the same as for Eq. (21). That is, the expression for the PDF given byEq. (42) is valid when ν > η > ν < η <
1. The average of a power of x , calculated using Eq. (42), is h x a i x = Z ∞ x a P ( x, t | x , dy = Γ (cid:16) ν − − a η − (cid:17) Γ (cid:16) ν − η − (cid:17) x a m (1 − e − η − µt ) a η − F a η −
1) ; ν − η −
1) ; − x η − x − η )0 e η − µt − ! . (43)This average is finite under the same conditions as Eq. (22). In particular, the average of x is equal to h x i x = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) x m (1 − e − η − µt ) η − F η −
1) ; ν − η −
1) ; − x η − x − η )0 e η − µt − ! (44)and is finite when ν > η > ν < η <
1. The average of the square of x is equal to h x i x = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) x (1 − e − η − µt ) η − F η −
1) ; ν − η −
1) ; − x η − x − η )0 e η − µt − ! (45)and is finite when ν > η > ν < η < µ has the same sign as η − t → ∞ then the PDF (42) tends to the steady-state PDF P ( x ) = 2 | η − | x ν − Γ (cid:16) ν − η − (cid:17) x − ν exp (cid:18) − (cid:16) x m x (cid:17) η − (cid:19) . (46)The steady-state PDF (46) leads to the steady-state averages of x and x h x i st = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) x m , (47) h x i st = Γ (cid:16) ν − η − (cid:17) Γ (cid:16) ν − η − (cid:17) x . (48)Now let us consider the time evolution of the average h x i x , given by Eq. (45). In the case when the initial position x is far from the cut-off boundary x m (that is x ≪ x m when η < x ≫ x m when η > h x i x can be separated into three parts. First, for small times t ≪ x − η )0 η − σ the influence of the initial position is significant and the diffusion is approximately normal, h x i x depends linearlyon time t . For the intermediate times x − η )0 η − σ ≪ t ≪ η − µ = x − η )m η − σ the exponent e − η − µt in Eq. (45) differs from 1 only slightly, however the last argument of the hypergeometricfunction is already small. Approximating the hypergeometric function by 1 and expanding the exponent e − η − µt into power series and keeping only the linear term we obtain that the average h x i x depends on time as a power-law, h x i x ∼ t / (1 − η ) . Thus for this intermediate range of time the anomalous diffusion occurs. Finally, for large times t ? η − µ the cut-off position x m starts to influence the diffusion and h x i x approaches the steady-state value (48). We canconclude that the introduction of the boundary via an exponential cut-off does not change the anomalous diffusionwhen the starting position is far from the boundary ant the time is not too large.Comparison of the analytic expressions (44), (45) with numerically obtained time-dependent mean and varianceis shown in Fig. 3. As in Sec. III, for numerical solution we use the Euler-Maruyama scheme with a variable timestep, equivalent to the introduction of the operational time in addition to the physical time t . When the diffusioncoefficient in the operational time τ does not depend on position x , the numerical method is given by the equations x k +1 = x k + (cid:20)(cid:16) η − ν (cid:17) x k + ( η − x η − x − ηk (cid:21) ∆ τ + √ ∆ τ ε k , (49) t k +1 = t k + ∆ τσ x ηk . (50)When η >
1, a more efficient numerical method is obtained when the change of the variable x in one step is proportionalto the value of the variable. Then the numerical method of solution is described by the equations x k +1 = x k " τ (cid:16) η − ν (cid:17) + ∆ τ ( η − (cid:18) x m x k (cid:19) η − + √ ∆ τ ε k , (51) t k +1 = t k + ∆ τσ x η − k . (52)0 -1 -1 ⟨ x ⟩ t (a) -1 -1 ⟨ x ⟩ t (b) -3 -2 -1 -1 ⟨ x ⟩ t (c) -2 -1 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (d) -2 -1 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (e) -4 -2 -1 ⟨ ( x - ⟨ x ⟩ ) ⟩ t (f) FIG. 3. Dependence of the mean (a,b,c) and variance (d,e,f) on time for various values of the parameters η and ν when theposition of the diffusing particle changes according to Eq. (40). Solid gray lines show numerical result, dashed black lines arecalculated using Eqs. (44) and (45), dotted (red) lines show the power-law dependence on time ∼ t / [2(1 − η )] for (a,b,c) and ∼ t / (1 − η ) for (d,e,f). The parameters are σ = 1 and η = − , ν = − x m = 5 for (a,d); η = , ν = 0, x m = 100 for (b,c); η = , ν = 5 x m = 0 .
01 for (c,f). The initial position is x = 1. -5 -4 -3 -2 -1 -2 -1 P ( x ) x (a) -5 -4 -3 -2 -1 -2 -1 P ( x ) x (b) -4 -2 -3 -2 -1 P ( x ) x (c) FIG. 4. Time-dependent PDF P ( x, t | x ,
0) corresponding to times t = 1 (light gray) and t = 10 (dark gray) for various valuesof the parameters η and ν when the position of the diffusing particle changes according to Eq. (40). Dashed black lines arecalculated using Eq. (42). The dotted line shows the steady-state PDF (46). The parameters are σ = 1 and (a) η = − , ν = − x m = 5; (b) η = , ν = 0, x m = 100; (c) η = , ν = 5 x m = 0 .
01. The initial position is x = 1. When η <
1, to avoid the divergence of the diffusion and drift coefficients at x = 0 in the numerical simulation, weinsert a reflective boundary at x = 10 − . This is not necessary when η > x ∼ x m . For numerical simulation we used the same values of theparameters η and ν as in Fig. 1, the initial position is x = 1. We have chosen the parameter x m of the externalforce so that the initial position x is far from the boundary x m . In Fig. 3 we see a good agreement of the numericalresults with analytic expressions. The numerical calculation and the analytic expressions confirm the presence of atime interval where the mean and the variance have a power-law dependence on time, as can be seen in Fig. 3. Theupper limit of this time interval is determined by the the position x m of the exponential cut-off.Comparison of the analytic expression (42) for the time-dependent PDF P ( x, t | x ,
0) with the results of numericalsimulation is shown in Fig. 4. We see a good agreement of the numerical results with the analytic expression. Withincreasing time the PDF shifts to the larger values of x when η < x when η > V. CONCLUSIONS
In summary, we have obtained analytic expressions (21) and (42) for the transition probability of the heterogeneousdiffusion process whose diffusivity has a power-law dependence on the distance x . In the description of the HDP wehave included an additional deterministic force that also has a power-law dependence on the position. The drift termhaving a power-law dependence on the position can arise not only due to an external force but can also representa noise-induced drift [45]. Such a drift term appears in a Langevin equation describing overdamped fluctuations ofthe position of a particle in nonhomogeneous medium [68]. Stochastic differential equations with power-law drift anddiffusion terms have been used to model random fluctuations of the atmospheric forcing on the ocean circulation [69]and pressure time series routinely used to define the index characterizing the North Atlantic Oscillation [70]. TheBrownian motion of a colloidal particle in water subjected to the gravitational force and with a space-dependentdiffusivity due to the presence of the bottom wall of the sample cell has been investigated in Ref. [71]. Force causingexponential cut-off in the PDF of the particle position can describe HDP process in confined regions; such a descriptioncan be relevant, for example, for the tracer particles moving in the confinement of cellular compartments [12, 20].A system obeying SDEs with power-law drift and diffusion terms can be experimetally realized as an electricalcircuit driven by a multiplicative noise, similarly as in Ref. [72]. The equation describing the overdamped motionof the Brownian particle takes the form of Eq. (10) when a temperature gradient is present in the medium and theparticle is subjected to the external potential that is proportional to the temperature profile [33]. In particular, steadystate heat transfer due to the temperature difference between the beginning and the end of the system corresponds to η = 1 / η −
1, where 2 η is the power-law exponent inthe dependence of the diffusion coefficient on the position, the external force does not limit the region of diffusion.Other values of the power-law exponent in the deterministic force can cause the exponential cut-off in the PDF of theparticle positions. Such an exponential restriction of the diffusion appears when the external force is positive if thepower-law exponent is smaller than 2 η − η −
1. We obtainedan analytic expression (42) for the transition probability in a particular case when the external restricting force hasa linear dependence on the position. Using analytic expression for the transition probability we calculated the timedependence of the moments of the particle position, Eqs. (22) and (43).We found that the power-law exponent in the dependence of the MSD on time does not depend on the externalforce, this force changes only the anomalous diffusion coefficient. This conclusion is valid for sufficiently large timessatisfying the condition (25), that is when the initial position of the particle is forgotten and the anomalous diffusionoccurs. In addition, the external force having the power-law exponent different from 2 η − α in Eq. (1) is determined only by the diffusioncoefficient.In addition, the results of Sec. III show that the character of the anomalous diffusion does not depend on theinterpretation of the Langevin equation: the scaling exponent α in Eq. (1) is the same for both Stratonovich and Itˆoconventions. This is because different interpretations correspond to the different values of the parameter ν in Eq. (10):Stratonovich convention results in ν = η , Itˆo convention in ν = 2 η , and the scaling exponent α does not depend on ν . The same conclusion that the exponent of the anomalous diffusion does not depend on the prescription has beenobtained in Refs. [58, 73, 74] for equations describing diffusion without the presence of an external force. Appendix A: Solution of the Fokker-Planck equation for the Bessel process
Let us consider the Fokker-Planck equation ∂∂t P = ν ∂∂x x P + 12 ∂ ∂x P (A1)and search for the time-dependent solution P ( x, t | x ,
0) with the initial condition P ( x, | x ,
0) = δ ( x − x ). Theunnormalized time-independent solution of Eq. (A1) is x − ν . The boundary condition at x = 0 for Eq. (A1) can beexpressed using the probability current [75] S ( x, t ) = − ν x P ( x, t ) − ∂∂x P ( x, t ) . (A2)2We consider the boundary condition corresponding to the vanishing probability current at x = 0, S (0 , t ) = 0.One of the possible ways to obtain the solution of Eq. (A1) is to use a Laplace transform [61]. Here we solveEq. (A1) using the method of eigenfunctions. This method has been used in Ref. [76] for an equation, similar toEq. (A1). An ansatz of the form P ( x, t ) = P λ ( x ) e − λt (A3)leads to the equation ν ∂∂x x P λ + 12 ∂ ∂x P λ = − λP λ , (A4)where P λ ( x ) are the eigenfunctions and λ > Z ∞ x ν P λ ( x ) P λ ′ ( x ) dx = δ ( λ − λ ′ ) . (A5)Expansion of the transition probability density P ( x, t | x ,
0) in terms of the eigenfunctions has the form [75] P ( x, t | x ,
0) = Z ∞ P λ ( x ) x ν P λ ( x ) e − λt dλ . (A6)For the solution of Eq. (A4) it is convenient to write the eigenfunctions P λ ( x ) as P λ ( x ) = x − γ u λ ( x ) , (A7)where γ = 1 + ν . (A8)Similar anzatz has been used in Refs. [62, 63]. The functions u λ ( x ) obey the equation x d dx u λ + x ddx u λ + ( ρ x − γ ) u λ = 0 , (A9)where ρ = √ λ . (A10)The probability current S λ ( x ), Eq. (A2), rewritten in terms of functions u λ , becomes S λ ( x ) = − x − γ (cid:18) γu λ ( x ) + x ddx u λ ( x ) (cid:19) . (A11)The orthonormality of eigenfunctions (A5) yields the orthonormality for functions u λ ( x ) Z ∞ xu λ ( x ) u λ ′ ( x ) dx = δ ( λ − λ ′ ) . (A12)The general solution of Eq. (A9) is u λ ( x ) = c J γ ( ρx ) + c Y γ ( ρx ) , (A13)where J γ ( x ) and Y γ ( x ) are the Bessel functions of the first and second kind, respectively. The coefficients c and c needs to be determined from the boundary and normalization conditions for the functions u λ ( x ). Using Eqs. (A13)and (A11) we get the probability current S λ ( x ) = − ρx − γ [ c J γ − ( ρx ) + c Y γ − ( ρx )] . (A14)3The requirement S λ (0) = 0 leads to the condition c = − c tan( πγ ). Taking into account this relation between c and c we obtain u λ ( x ) = c λ J − γ ( ρx ) . (A15)In addition, the condition S λ (0) = 0 implies c = 0 when γ >
1. Thus, when γ > ν > c and c are zero and the solution (A6) is not valid. It is known that for a Bessel process with such parameters a totalabsorption at the origin occurs in a finite time [64].Orthonormality condition (A12) leads to the equation c λ c λ ′ Z ∞ xJ − γ ( ρx ) J − γ ( ρ ′ x ) dx = δ ( λ − λ ′ ) . (A16)Since for the Bessel functions the equality ρ Z ∞ xJ γ ( ρx ) J γ ( ρ ′ x ) xdx = δ ( ρ − ρ ′ ) (A17)is valid, we obtain c λ = 1. Using Eqs. (A6), (A7) and (A15) the solution of the Fokker-Planck equation can beexpressed as P ( x, t | x ,
0) = x − γ x γ Z ∞ J − γ ( ρx ) J − γ ( ρx ) e − ρ t ρdρ . (A18)Integration yields P ( x, t | x ,
0) = x − γ x γ t exp (cid:18) − x + x t (cid:19) I − γ (cid:16) xx t (cid:17) . (A19)Here I α ( x ) is the modified Bessel function of the first kind. [1] J.-P. Bouchaud and A. Georges, Phys. Rep. , 127 (1990).[2] F. Evers, C. Zunke, R. D. L. Hanes, J. Bewerunge, I. Ladadwa, A. Heuer, and S. U. Egelhaaf, Phys. Rev. E , 022125(2013).[3] O. Chepizhko and F. Peruani, Phys. Rev. Lett. , 160604 (2013).[4] R. Metzler and J. Klafter, Phys. Rep. , 1 (2000).[5] M. Schubert, E. Preis, J. C. Blakesley, P. Pingel, U. Scherf, and D. Neher, Phys. Rev. B , 024203 (2013).[6] C. Scalliet, A. Gnoli, A. Puglisi, and A. Vulpiani, Phys. Rev. Lett , 198001 (2015).[7] H. C. Fogedby, Phys. Rev. Lett. , 2517 (1994).[8] Y. Marandet, H. Capes, L. Godbert-Mouret, R. Guirlet, M. Koubiti, and R. Stamm, Commun. Nonlinear. Sci. Commun. , 469 (2003).[9] N. Mercadier, W. Guerin, M. Chevrollier, and R. Kaiser, Nat. Phys. , 602 (2009).[10] S. M. A. Tabei, S. Burov, H. Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L. H. Philipson, A. R. Dinner, and N. F.Scherer, Proc. Natl Acad. Sci. USA , 4911 (2013).[11] A. V. Weigel, B. Simon, M. M. Tamkun, and D. Krapf, Proc. Natl Acad. Sci. USA , 6438 (2011).[12] J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys.Rev. Lett. , 048103 (2011).[13] I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch, and D. A. Weitz, Phys. Rev.Lett. , 178101 (2004).[14] H. Scher and E. W. Montroll, Phys. Rev. B , 2455 (1975).[15] J.-H. Jeon, H. Martinez-Seara Monne, M. Javanainen, and R. Metzler, Phys. Rev. Lett. , 188103 (2012).[16] G. R. Kneller, K. Baczynski, and M. Pasienkewicz-Gierula, J. Chem. Phys , 141105 (2011).[17] E. Kepten, I. Bronshtein, and Y. Garini, Phys. Rev. E , 052713 (2013).[18] J. Szymanski and M. Weiss, Phys. Rev. Lett. , 038102 (2009).[19] J.-H. Jeon, N. Leijnse, L. B. Oddershede, and R. Metzler, New J. Phys. , 045011 (2013).[20] T. K¨uhn, T. O. Ihalainen, J. Hyv¨aluoma, N. Dross, S. F. Willman, J. Langowski, M. Vihinen-Ranta, and J. Timonen,PLoS ONE , e22962 (2011).[21] A. G. Cherstvy, A. V. Chechkin, and R. Metzler, New J. Phys. , 083039 (2013).[22] A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys , 20220 (2013).[23] A. G. Cherstvy, A. V. Chechkin, and R. Metzler, Soft Matter , 1591 (2014). [24] A. G. Cherstvy and R. Metzler, Phys. Rev. E , 012134 (2014).[25] Y. T. Maeda, T. Tlusty, and A. Libchaber, Proc. Natl. Acad. Sci. , 17972 (2012).[26] C. B. Mast, S. Schink, U. Gerland, and D. Braun, Proc. Natl. Acad. Sci. , 8030 (2013).[27] R. Haggerty and S. M. Gorelick, Water Resources Res. , 2383 (1995).[28] M. Dentz and D. Bolster, Phys. Rev. Lett. , 244301 (2010).[29] L. F. Richardson, Proc. R. Soc. Lond. A , 709 (1926).[30] B. O’Shaughnessy and I. Procaccia, Phys. Rev. Lett. , 455 (1985).[31] C. Loverdo, O. B´enichou, R. Voituriez, A. Biebricher, I. Bonnet, and P. Desbiolles, Phys. Rev. Lett. , 188101 (2009).[32] B. P. English, V. Hauryliuk, A. Sanamrad, S. Tankov, N. H. Dekker, and J. Elf, Proc. Natl Acad. Sci. USA , E365(2011).[33] R. Kazakeviˇcius and J. Ruseckas, J. Stat. Mech. , P02021 (2015).[34] T. Srokowski, Phys. Rev. E , 030102(R) (2014).[35] D. Schertzer, M. Larchevˆeque, J. Duan, V. V. Yanovsky, and S. Lovejoy, J. Math. Phys. , 200 (2001).[36] M. Dentz, P. Gouze, A. Russian, J. Dweik, and F. Delay, Adv. Water Res. , 13 (2012).[37] A. E. Cohen, Phys. Rev. Lett. , 118102 (2005).[38] G. S. Manning, J. Chem. Phys. , 924 (1969).[39] Y. Sagi, M. Brook, I. Almog, and N. Davidson, Phys. Rev. Lett. , 093002 (2012).[40] F. Bouchet and T. Dauxois, Phys. Rev. E , 045103(R) (2005).[41] L. A. Wu, S. S. Wu, and D. Segal, Phys. Rev. E , 061901 (2009).[42] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Springer-Verlag, Berlin,2004).[43] A. W. C. Lau and T. C. Lubensky, Phys. Rev. E , 011123 (2007).[44] A. Fuli´nski, Phys. Rev. E , 061140 (2011).[45] G. Volpe and J. Wehr, Rep. Prog. Phys. , 053901 (2016).[46] B. Kaulakys and J. Ruseckas, Phys. Rev. E , 020101(R) (2004).[47] B. Kaulakys, J. Ruseckas, V. Gontis, and M. Alaburda, Physica A , 217 (2006).[48] V. Gontis, J. Ruseckas, and A. Kononovicius, Physica A , 100 (2010).[49] J. Mathiesen, L. Angheluta, P. T. H. Ahlgren, and M. H. Jensen, Proc. Natl. Acad. Sci. , 17259 (2013).[50] R. Ton and A. Daffertshofer, NeuroImage (2016), doi: 10.1016/j.neuroimage.2016.01.008.[51] I. Karatzas and S. Shreve, Brownian Motion and Stochastic Calculus (Springer, New York, 2012).[52] P. H¨anggi and H. Thomas, Phys. Rep. , 207 (1982).[53] Y. L. Klimontovich, Phys. Usp. , 737 (1994).[54] G. Pesce, A. McDaniel, S. Hottovy, J. Wehr, and G. Volpe, Nat. Commun. , 2733 (2013).[55] B. C. dos Santos and C. Tsallis, Phys. Rev. E , 061119 (2010).[56] I. M. Sokolov, Chem. Phys. , 359 (2010).[57] N. G. van Kampen, J. Stat. Phys. , 175 (1981).[58] M. Heidern¨atsch, On the diffusion in inhomogeneous systems , Ph.D. thesis, Technische Universit¨at Chemnitz, Faculty ofSciences, Institute of Physics, Complex Systems and Nonlinear Dynamics (2015).[59] S. F. Kwok, Ann. Phys. , 1989 (2012).[60] Z. G. Arenas, D. G. Barci, and C. Tsallis, Phys. Rev. E , 032118 (2014).[61] M. Jeanblanc, M. Yor, and M. Chesney, Mathematical Methods for Financial Markets (Springer, London, 2009).[62] A. J. Bray, Phys. Rev. E , 103 (2000).[63] E. Martin, U. Behn, and G. Germano, Phys. Rev. E , 051115 (2011).[64] S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes , 2nd ed. (Academic, New York, 1981).[65] J. Ruseckas, R. Kazakevicius, and B. Kaulakys, J. Stat. Mech. , 054022 (2016).[66] A. G. Cherstvy, A. V. Chechkin, and R. Metzler, J. Phys. A: Math. Theor. , 485002 (2014).[67] A. G. Cherstvy and R. Metzler, J. Stat. Mech. , P05010 (2015).[68] J. M. Sancho, M. San Miguel, and D. D¨urr, J. Stat. Phys. , 291 (1982).[69] P. D. Ditlevsen, Geophys. Res. Lett. , 1441 (1999).[70] P. G. Lind, A. Mora, J. A. C. Gallas, and M. Haase, Phys. Rev. E , 056706 (2005).[71] G. Volpe, L. Helden, T. Brettschneider, J. Wehr, and C. Bechinger, Phys. Rev. Lett. , 170602 (2010).[72] J. Smythe, F. Moss, and P. V. E. McClintock, Phys. Rev. Lett. , 1062 (1983).[73] K. S. Fa and E. K. Lenzi, Phys. Rev. E , 061105 (2003).[74] K. S. Fa, Phys. Rev. E , 020101 (2005).[75] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer-Verlag, Berlin, 1989).[76] J. Ruseckas and B. Kaulakys, Phys. Rev. E81