Monte Carlo evaluation of FADE approach to anomalous kinetics
1 Monte Carlo Evaluation of FADE Approach to Anomalous Kinetics
M. Marseguerra , A. Zoia Department of Nuclear Engineering, Polytechnic of Milan, Via Ponzio 34/3, 20133 Milan, Italy
Abstract
In a wide range of transport phenomena in complex systems, the mean squared displacement of a particles plume has been often found to follow a nonlinear relationship of the kind a ttx (cid:181) )( , where a may be greater or smaller than : these evidences have been described under the generic term of anomalous diffusion. In this paper we focus on subdiffusion, i.e. the case , in presence of an external advective field. Widely adopted models to describe anomalous kinetics are Continuous Time Random Walk (CTRW) and its Fractional Advection-Dispersion Equation (FADE) asymptotic approximation, which accurately account for experimental results, e.g. in the transport of contaminant particles in porous or fractured media. FADE approximated equations, in particular, admit elegant analytical closed-form solutions for the particle concentration P(x,t) . To evaluate the relevance of the approximations which allow to derive the asymptotic FADE equations, we resort to Monte Carlo simulation (which may be regarded as an exact solution of the CTRW model): this comparison shows that the FADE equations represent a less and less accurate asymptotic description of the exact CTRW model as a becomes close to . We propose higher order corrections which lead to modified integral-differential equations and derive new expressions for the moments of P(x,t) . These results are validated through comparison with those of Monte Carlo simulation, assumed as reference curves.
1. Introduction
A general approach to the analysis of transport phenomena is based on Continuous-Time Random Walk (CTRW) [15,19,22,23,35], in which the travel of a particle (a walker) in a medium is modelled as a series of jumps of random lengths, separated by random waiting times. The theory of CTRW with algebraically decaying probability distribution functions (pdf’s) has been originally introduced in Physics in a series of seminal papers by Weiss, Scher, Montroll and co-workers [15,27,29,30,31] in the late 1960s to explain evidences of anomalous diffusion occurring in the drift-diffusion processes in amorphous semiconductors. The diffusion is called anomalous if the mean squared displacement (MSD) of the particles plume is not linearly proportional to time t as in the standard Fickian case, but to a power a of t larger (superdiffusion) or smaller (subdiffusion) than unity. More recently, anomalous diffusion has turned out to be quite ubiquitous in almost every field of science (see e.g. [22] and [23] for a detailed review) and the CTRW model has been applied with success to interpret the experimental results and to make predictions on the evolution of the examined systems. Such applications concern among the others e.g. the behaviour of chaotic Hamiltonian systems as related to the transport of charged particles in turbulent plasma [33,36], the evolution of financial markets [18], the dynamics of disordered systems [7] or the transport of contaminant particles in groundwater under the combined effect of rock fractures and porosity [3-6,8-10,16,20,21,25]. In this paper we will focus on the case of subdiffusion ( ) on a 1D infinite support, in presence of an external (constant) advective field. In the subdiffusive case [3-6,8-10,16,19-21,22,25], an algebraically decaying distribution is assumed for the waiting times of the particles in the surrounding medium (instead of the traditional exponential one): this physically means that the particles will have a non-vanishing probability of extremely long sojourn times in the visited locations, due to so-called “trapping events” (see e.g. [36]). The macroscopic effect is that the variance of the particle plume grows less than linearly in time (subdiffusion). Within the CTRW approach to transport phenomena, advection can be modelled in either a Galilei invariant [22,25] or Galilei variant [3,4,20,22,25] scheme, each involving different hypotheses on the microscopic dynamics of the described particle plume. In general, the CTRW equations do not allow closed-form analytical solutions. However, if suitable first-order approximations are introduced in the Fourier and Laplace transforms expansions for both the jump lengths Corresponding Author. Tel: +39 02 2399 6355. Fax: +39 02 2399 6309 Email address: [email protected] (Marzio Marseguerra) Email address: [email protected] (Andrea Zoia) and waiting times distributions, Fractional Advection-Diffusion Equations (FADE) can be formally derived from the original exact CTRW model [22,23,32,33,34,36]. In this respect, FADE may be regarded as an asymptotic subset of CTRW. In order to evaluate the effects of the approximations introduced in the transformed space as to obtain the FADE equations, we resort to Monte Carlo simulation as a mean to accurately solve CTRW equation (as proposed e.g. in [3,10,21]): it turns out that relevant deviations of the analytical FADE results with respect to the Monte Carlo solutions are evident when fia . In particular, our analysis focuses on the moments of the distribution P(x,t) , which constitute accurate estimators of the solution accuracy [19,21]. These observations suggest that when fia FADE is not an accurate asymptotic solution of CTRW and show the need of improving the first-order expansions which lead to the FADE equations in order to include higher-order terms: in this case, no analytical solutions are available and numerical Laplace inversion will be required to obtain the desired solutions of the CTRW equation and their moments. The proposed corrections to the transformed space expansions, validated by means of Monte Carlo simulation, show that the new set of derived equations can accurately approximate the behaviour of the exact CTRW solution even when fia . The paper is organized as follows: in Section 2 we summarize the basic concepts of the CTRW approach and the derivation of the Fractional Diffusion Equation (FDE). In Section 3 we present the Galilei invariant and variant schemes within the CTRW approach and we derive the corresponding FADE: analytical results of the asymptotic FADE will be compared to Monte Carlo simulation results. We will show that discrepancies arise when fia . In Section 4 we propose higher order corrections for both advection schemes and compare the numerical results to the Monte Carlo simulation findings. Conclusions are finally drawn in Section 5. Appendices A and B are devoted to the required calculations.
2. CTRW and FDE approaches to subdiffusion
Let
X(t) be a stochastic process describing the motion of a tracer particle (a walker) performing random jumps separated by random waiting times. The associated probability density function (pdf)
P(x,t) represents the probability of the walker being in
X=x at time t and it is also called the propagator of the process. The CTRW approach to the description of this stochastic process is based on the probability balance expressed by the Chapman-Kolmogorov integral equation (the so-called Master Equation), which entails the pdf P(x,t) [22,35,36]. It can be shown [12,18,22,23,36] that, if )( x l and )( tw are the pdf’s of the single jump lengths and waiting times , respectively, the Laplace- and Fourier- transformed expression of P(x,t) satisfies the simple algebraic relation )()(1 1)(1),( kuwu uwukP l--= (1) in the case of a Cauchy problem with initial conditions )()( tx dd . For convenience, adopting a well-established convention (see e.g. [19,22,23]), we denote the Fourier or Laplace transform of a pdf by its argument, namely { } )();( kfkxf =` and { } )();( ufutf = L . P(k,u) is called the propagator of the underlying stochastic process: equation (1), in the doubly transformed space, represents a probability (normalized mass) balance for the number of particles and in this respect may be regarded as an exact transport formulation. According to the choice of the waiting times and of the jump lengths distributions, the CTRW approach may give rise to (normal) diffusive, subdiffusive and superdiffusive behaviour of the walker. For the case of subdiffusion, i.e. , the most widely adopted choice for )( x l is a Gaussian pdf with finite variance, while )( tw is assumed to be any algebraically decaying pdf of the kind a-- ~)( ttw when +¥fi t . The specific functional form of )( tw , whose first moment is infinite, introduces “memory effects” in the trajectory )( txx = of the walker, giving rise to long range correlations which make its path semi-Markovian . This fact in turn prevents the application of the Central Limit Theorem and results in anomalous diffusion, i.e. a ttx (cid:181) )( [11,18,22,36]. In general, as mentioned before, no analytical solution is known for (1), as the required Laplace and Fourier inverse transforms are usually not trivial . However, if we consider the asymptotic behavior of the solution sufficiently far from the origin, i.e. the so-called “diffusion limit” approximation { } { } fifi«+¥fi+¥fi uktx and expand the Laplace and Fourier transforms )( k l and Here, for sake of simplicity, we adopt the hypothesis that the jump length and waiting time probabilities are independent. In general, however, they may well be correlated. See below and Appendix B for details on the semi-Markovian nature of the walker when a<1 . It must be mentioned here that it is possible to resort to a rephrasing of the original probability balance Master Equation in terms of an ordinary differential equation which can be then analytically solved with respect to space to get P(x,u). However, to obtain the full solution of the CTRW equation P(x,u) needs to be numerically Laplace inverted with respect to time [8-10]. )( uw up to the first non-constant term, the required inverse transforms may be explicitly performed. Therefore, from (1) by means of the definition of the differential operators in the transformed spaces we can formally derive a Fractional Diffusion Equation (FDE) whose analytical solution is known and may be expressed in terms of the Fox H function, as first shown by Schneider and Wyss [17,22,26,32]. In the following, for sake of simplicity we will assume a pdf for the waiting times of the kind taat tt >+£ += tt tAtAtw )( (2) where A is a normalization factor, which reads aa+= A , and t is a characteristic time constant. This pdf has an asymptotic behaviour a+ » ttw when +¥fi t , so that we expect w(t) to converge to a one-sided Lévy stable distribution of index a [11]. Indeed, if the Laplace transform of (2) is expanded in the long time limit (i.e. < 3. Modelling an advection field within CTRW approach Once the general form of the propagator (1) has been set, there basically exist two schemes which allow to take into account an external velocity field in presence of subdiffusion: this distinction has been first introduced in [25]. Here, for sake of simplicity, we will consider a constant homogeneous advective field. The first scheme is introduced by requiring that the solution P(x,t) of (1) be invariant under a transformation of coordinates of the kind vtxx -fi : this idea stems directly from the standard Fickian Advection-Dispersion Equation (ADE). This scheme, which is called “Galilei invariant” because of the mentioned invariance property, assumes that in the moving frame (at velocity v ) each particle is ruled by the usual pdf’s )( x l and )( tw . Galilei invariant subdiffusion may be found e.g. for tracer particles immersed in polymer solutions, where the flowing substance itself is the cause of the slowing down of the motion [22,25]. If we indicate ),( tx V and ),( tx c the jump probability in the moving frame and in the laboratory frame, respectively, then the assumption of Galilei invariance implies that ),(),( tvtxtx -= Vc , or equivalently ),(),( ivkukuk += Vc [22,25]. (6) If we assume expansions (3) and (4) to hold true, propagator (1) readily becomes kuDivkuukP aa - ++= (7) Again, recalling the properties of a- ¶ t in the transformed space and with the help of some algebra, propagator (7) may be formally inverted to the { } tx , space to give ),(),(),( txPxDtxPxvtxPt t ¶ ¶¶=¶ ¶+¶ ¶ - aa (8) which is called the (Galilei invariant) Fractional Advection Dispersion Equation (FADE) and may be considered as a generalization of the well-known ADE. Its analytical solution may be obtained in terms the solution of (5) by observing that by definition )0;,();,( =-= vtvtxPvtxP . Moreover, recalling that by definition ( ) ( ) ¶ ¶= fi ),(lim ukPkitx nnnk-n L , (9) from equation (7) it is possible to deduce all the moments of P(x,t) : in particular, under approximations (3) and (4), the Laplace inverse transform in (9) can be analytically computed, to obtain ( ) vttx = , (10) ( ) ( ) )1(2 vttDtx ++G= aa a , (11) so that ( ) ( ) aa a tDtxtx )1(2)( +G=- . (12) The moments of distribution P(x,t) are sensitive estimators of the accuracy of P(x,t) in (8) as an approximation of the exact solution of (1). Since, as required by the Galilei invariance property, the variance (12) has exactly the same expression of the one of the FDE without advection [22,25], the accuracy of the FADE (8) is traced back to the one of the FDE and we expect expressions (11) and (12) to be accurate when a <0.5 and to show relevant discrepancies as fia . The first moment (10) is not influenced by the order of expansion (3). In the following Figures for two different values of a we compare the analytical approximate variance (12) with the one computed by resorting to Monte Carlo simulation as described in [21]. The simulation parameters are the following: 10 particles have been followed up to a time t max = 10 , with s =0.5, t =1 and v =1 when a =0.5 (Figure 1). Then we set s =0.5, t =1, v =1 and t max =10 when a =0.8 (Figure 2). We remark that in both cases the explored time scales guarantee that t is “small”, thus ensuring that the diffusion limit is a justified assumption. As expected, for large a FADE (8) is not an accurate approximation of (1) (Figure 2). Figure 1. Comparison between analytical variance (12) (solid line) and Monte Carlo variance (dots) when a =0.5. Figure 2. Comparison between analytical variance (12) (solid line) and Monte Carlo variance (dots) when a =0.8. We proceed now further to examine the other possible scheme to model an advection field within the CTRW approach, namely the Galilei variant scheme. Starting again from (1), it is assumed that each jump has a spatial (constant) bias in the forward direction, m . This assumption does not result in a Galilei invariance property for P(x,t) , hence the name of the scheme [25]. Experimental evidences of this kind of behaviour may be e.g. found in the context of flow through porous or fractured media (and in general in disordered systems) [3,4,19,20,22,25]. This scheme is modelled by adding to the jump lengths pdf )( x l a finite mean m . Thus, the expansion in small k of the Fourier transform (4) becomes kkik sml -+» . (13) We remark that m has the dimensions of a space. Assuming that expansion (3) still holds true, one readily deduces from (1) the propagator aaaa -- ++= kuivukDuukP , (14) where aaa tm cv = is a generalized advection coefficient (or velocity). Again, by taking into account the properties of the Riemann-Liouville operator, propagator (14) may be formally inverted in the { } tx , space, to finally obtain ¶ ¶-¶ ¶¶=¶ ¶ - ),(),(),( txPxvtxPxDtxPt t aaa . (15) Equation (15) is called (Galilei variant) FADE or Fractional Fokker-Plank-Kolmogorov Equation (FFPK) [18,33,36]: its structure is neatly different from that of the Galilei invariant FADE, because of the presence of the fractional differential operator acting also on the advective contribution. Now, by making use of definition (9), we can derive all the moments of the distribution P(x,t) solution of (15): their comparison with those computed by means of Monte Carlo simulation will provide an indication on the accuracy of the asymptotic P(x,t) with respect to the exact solution of CTRW (1). After some algebra, we obtain: ( ) )1( a aa +G= tvtx (16) ( ) ( ) aaaa aa )21(2)1(2 tvtDtx +G++G= . (17) In the following figures for two different values of a we compare the analytical approximate variance ( ) ( ) )( txtx - with the one computed by resorting to Monte Carlo simulation as described in [21]. The simulation parameters are the following: 10 particles have been followed up to a time t max = 10 , with s =0.5, t =1 and m =0.2 when a =0.5 (Figure 3). Then we set s =0.5, t =1, m =0.2 and t max =10 when a =0.8 (Figure 4). As expected, for large a FADE (15) is not an accurate approximation of (1) (Figure 4). Figure 3 Comparison between analytical (solid line, from (16) and (17)) and Monte Carlo variance (dots) when a =0.5. Figure 4. Comparison between analytical (solid line, from (16) and (17)) and Monte Carlo variance (dots) when a =0.8. 4. Higher order expansion of Laplace transform (3) In reason of the considerations presented in the previous paragraphs and elsewhere [19,21], we can attribute the discrepancies shown in Figures 2 and 4 to the truncation to the first non-constant term in (3), which in turn have led to the FADE formulation. Indeed, it can be shown that the Laplace transform expansion of pdf (2) in the diffusion limit < 112 212 cucutxtx - L (21) for the variance of (19) and -= - tt m aaa cucutx - L (22) ( ) -+ -= - ttstt m aaaaaa 112 2121212 cucuucucutx -- LL (23) for the first and second moments of (20). The inverse Laplace transform appearing in (21-23) may be computed with a numerical algorithm [13-14] and the obtained curves are compared with those obtained from the Monte Carlo simulations described above: the next Figures 5 and 6 show the improved accuracy of the second order asymptotic equations (21-23) with respect to expressions (12,16-17). The accuracy of the proposed higher order corrections of expansion (3) are fully validated by the comparison with the Monte Carlo results, which may be regarded as a reference solution for the CTRW, no approximation having been introduced in the stochastic modelling, but for finite statistics effects. Figure 5. Comparison between numerically inverted variance (21) and Monte Carlo variance (dots) when a =0.8. Figure 6. Comparison between numerically inverted variance (from (22) and (23)) and Monte Carlo variance (dots) when a =0.8. 5. Conclusions In a wide range of complex systems, the relaxation following an initial pulse is experimentally found to follow a nonlinear relationship of the kind a ttx (cid:181) )( , where a may be greater or smaller than : these evidences have been described under the generic term of anomalous kinetics. In particular, in this paper we have focused on the case of subdiffusion, i.e. In order to properly characterize the time evolution of such systems, several approaches have been proposed in literature, among which the CTRW model has been shown to be particularly successful. In particular, considering the so-called “diffusion limit” behavior of the evolving system (i.e. both spatially and temporally far from the source), the CTRW transport equation for a system undergoing anomalous diffusion in presence of an external advective field may be approximated by the asymptotic FADE model. This approximation is derived under the further assumption of truncating a required Laplace transform to the first non-constant term. The FADE model, by making use of the properties of fractional order derivatives, allows to obtain elegant close form analytical solutions. In this paper, by means of Monte Carlo simulation, we have explored the limits of validity of the FADE asymptotic solutions with respect to those of the CTRW, which must be in general obtained numerically. Both the case of a Galilei invariant and a Galilei variant advection schemes have been considered. Comparison with Monte Carlo results has revealed that the FADE solutions accurately represent the asymptotic behavior of the CTRW model only when a is far from its superior limit 1. As in most experimental evidences of anomalous kinetics a is usually close to 1 [see e.g. 8,18], we believe that the accuracy of the FADE solution should be improved. To this aim, we have expanded the cited Laplace transform by including higher order terms, whose importance has been shown to grow when a approaches 1. The proposed corrections do not allow in general to derive closed form solutions anymore, so that numerical inversion has been necessary: the solutions thus obtained turned out to be in very good agreement with the Monte Carlo solutions for CTRW, which have been adopted as reference curves. An appendix has been devoted to the derivation of the integral-differential equations corresponding to the improved Laplace transform expansions. Appendix A: Derivation of the Laplace transform of (2) In order to obtain the Laplace transform of (2), we proceed as follows: by definition, ∫∫∫ +¥ ----¥ - +++== t aat taataa dttedttedttweuw tututu 10 20 . By changing the variable of integration sut = , we get ∫∫∫ +¥ ----¥ - +++== t aaat taataa u su stu dsseudsu sedttweuw 10 220 . Then, by integrating by parts, the first integral reads ( ) ( ) ttaataa tt ueudsu se uu s +-+=+ -- ∫ 220 22 . Expanding the exponential in the diffusion limit <