Probability distribution functions of sub- and super-diffusive systems

We study the anomalous transport in systems of random walks (RW's) on comb-like lattices with fractal sidebranches, showing subdiffusion, and in a system of Brownian particles driven by a random shear along the x-direction, showing a superdiffusive behavior. In particular, we discuss whether scaling and universality are present or not in the shapes of the particle distribution along the preferential transport direction (x-axis).


I. INTRODUCTION
Standard diffusion is a Gaussian behavior characterized by a linear time growth of the mean square displacement (MSD) from the initial condition However certain processes in nature as well as in finance and even in sociology, are not Gaussian showing a nonlinear time growth of the MSD: If ν < 1/2 the process x(t), by definition, undergoes a sub-diffusive behavior while, if ν > 1/2, x(t) is said to be super-diffusive. There is not a unique origin of the anomalous diffusion; it can be due to many different causes which are specific to the process under investigation.
From the mathematical side, anomalous diffusion is a strict consequence of the breakdown of the Central Limit Theorem hypothesis, which can be ascribed either to the broad distributions of independent elementary steps, e.g. Lévy flights [1], or to the emergence of long-range, spatial or temporal, correlations among such steps.
In transport processes, the main subject of this paper, spatial correlations may arise from strong geometrical confinement and inhomogeneity due to the presence of disorder, obstacles, compartments and trapping sites, that can be found in amorphous [2,3] and porous [4,5] materials. Similarly, the crowding of cellular cytoplasm [6,7] and organic tissues [8,9] can make the motion of molecules and water in biological environments strongly correlated. Correlations among consecutive displacements can also emerge dynamically [10], as in the case of chaotic dynamics [11][12][13], or because the particles are driven by the motions of the underlying medium, e.g., diffusion in turbulent fluid. Richardson, for instance, proved that in a turbulent fluid, spatial dispersion of particles follows a super-diffusive behavior [14,15].
Similar situations occur when the random walks (RWs) are restricted on peculiar topological structures, such as fractals and networks [16][17][18], where the geometrical constraints do not allow a fast memory loss among a series of consecutive displacements.
Finally, it is important to remark that often, in real physical and biological cases [19], the anomalous diffusion can manifest itself as a transient phenomenon that, although long-lived, soon or later either turns into a standard one, or even stops due to the finiteness of the environments.
The goal of this paper is to discuss the general scaling properties of the probability distribution function (PDF) of anomalous diffusion at large times. In analogy with standard diffusion which undergoes the natural scaling we wonder to what extent the conjecture can hold for sub-and super-diffusive processes. Obviously, the consequent issue is to determine the shape of F ν (z) and understanding its degree of universality. At first, we have to bear in mind that the universality of this scaling is certainly broken by the existence of diffusion phenomena classified as "strongly anomalous" [13,20,21]. Indeed, they exhibit a multiscaling in the spectrum of moments, with not constant ν(m), which is in contrast with Eq.(2) because a single exponent is not sufficient to characterize the statistical properties of such processes. The conflict is evident by considering that Eq.(2) prescribes the following scaling for the moments, In the following we will not address strongly anomalous phenomena. Roughly speaking, there are two different scenarios that can be traced back to the property (2): one associated with time-homogeneous diffusion processes, for which the increment x(t+h)−x(t) is independent of x(t), and the other associated with processes not fulfilling the above property. The latter are said non-homogeneous processes and are especially observed in turbulent diffusion.
In the literature, some guesses can be found as to the functional form of F ν (z) depending on the anomalous exponent ν. It is worth mentioning two classes of F ν (z), the one suggested for the time-homogeneous processes (Fisher) and the one valid for turbulent diffusion (Richardson) both characterized by a stretched exponential with exponents Figure 1 sketches the expected validity range of the Fisher and Richardson F ν (z) as a function of the anomalous exponents, in both sub-and super-diffusive regimes. It is interesting to note that for ν = 1/2, the two stretched exponentials become of Gaussian type, as α = 2. We refer to the F ν (z) with α = 1/(1−ν) as the Fisherlike distribution, since Fisher derived it in the context of self-avoiding polymers (SAP) [23]. Whereas, the case of F ν (z) with α = 1/ν, which can be referred to as Richardson's function because it was derived by Richardson in his studies on turbulent diffusion [14,15], will not be addressed in this work. Before discussing the numerical analysis we carried out for proving the Fisher's scenario, we present two arguments supporting the validity of One is based on simple probability considerations [24], and the other on the large-deviation theory (LDT) which has been already and successfully applied to describe the continuous time RW (CTRW) [25,26].
The first heuristic argument starts from the observation that the main contribution to the tails arises from RWs that are very persistent along the x-direction, i.e. x max ∼ t. Therefore, the probability of such persistent walks can be approximated as On the other side, the probability of such paths is also P {z max } ∼ exp(−pt), as they undertake n = t/∆t independent steps in the same direction, p being a constant and ∆t the time step. Equating F ν (z max ) = P {z max } gives the Fisher's branch scaling prediction.
The second heuristic argument is based on LDT, which is a natural probability framework to determine how large fluctuations from the mean of a random process characterize the decay of far tails of its PDF [27]. LDT is the proper mathematical tool if we are interested in verifying whether the far tails of the PDF (2) are consistent with Fisher's tails. LDT assumes that, at large time, the probability decays as where C(· · · ) is the Cramers' (or rate) function. We recall that if a process µ n is the average over n independent random variables, µ n = (ξ 1 + · · · + ξ n )/n, its Cramers' function is, by definition, the limit C(µ n ) = − lim n→∞ 1/n ln P (µ n , n), implying that asymptotically for large n This expression can be equivalently rewritten in terms of the sum S = ξ 1 + · · · + ξ n = n µ n and the total time t = n∆t, which is exactly Eq. (8). A comparison between Eq. (8) and Eq.(5) implies that tC(x/t) = a ν |x/t ν | α , this equality can be rearranged to C(x/t) = a ν |x/t (ν+1/α) | α . Now the only consistent possibility is that C(x/t) = a ν |x/t| α which can only hold true if Eq.(7) is verified. Of course, although the above reasoning does not constitute a proof, yet we can conclude that Fisher's PDF is the unique PDF which is consistent with the scaling (2) and the LDT as well.
In this paper we will focus only on the validity of the Fisher's scenario in the subdiffusive and superdiffusive regimes, by using two models.
(1) The random walks (RWs) on a class of comb-like structures consisting of a main backbone decorated by an array of fractal sidebranches [28], as in Fig.2(A). In this case, the subdiffusion is observed along the backbone.
(2) The Lagrangian dynamics of particles in a channel geometry under the combined effects of a random velocity shear and molecular diffusivity, Fig.2(B) [22,29]. In this case, the superdiffusion occurs along the longitudinal transport direction (x-axis). We will show that the scaling behavior (2) is well satisfied by both the sub-and super-diffusive processes that we studied. This is testified by the perfect collapse, upon re-scaling x → |x|/t ν , of the simulated particle distribution (PDF) onto a master curve whose tails can be described by Eq.(5) with exponent (7). As a matter of fact we will demonstrate, numerically and analytically, that the Fisher relation is rigorously satisfied in the case of comb-like structures, even when the sidebranches are fractal structures of growing complexity (Sierpinski gaskets with increasing spectral dimension, Fig.2(A).
On the other side, the superdiffusive case of the random-shear model [ Fig.2(B)] presents an unexpected scenario where the Fisher scaling appears to be satisfied only for a restricted range of the anomalous exponent.
Moreover, the behavior of the PDF moments in agreement with Eq.(4) provides further numerical support to the validity of the re-scaling, Eq.(2).
As we will see in the following, the very numerical difficulty in observing the true nature of the tails stems from the necessity of two simultaneous asymptotic conditions, t 1 and also |x|/t ν 1. This paper is organized as follows. In Sec.II we study the anomalous behavior of RWs on comb-like structures with fractal sidebranches. In Sec.III, we analyze the numerical results on the super-diffusion of Brownian particles under a random shear in a channel of width L.
Finally Sec.IV contains a discussion and conclusions.

II. COMB-LIKE SYSTEMS
To study the scaling of the PDF of anomalous subdiffusion, we considered RW's on comb-like structures, namely, a central backbone decorated by either linear or more complex sidebranches (SB's) [17]. In particular, we are interested in fractal structures, so we consider SB's constituted of Sierpinski-gaskets with increasing spectral dimension [see Fig.2(A)]. Such branched topology is typical of percolation clusters at criticality, which can be viewed as finitely ramified fractals [30,31]. Comb-like structures, moreover, are frequently used in condensed matter and biological frameworks to describe the geometry of branched polymers [32,33], amphiphilic molecules and engineered structures at the nano-and microscales, and anomalous propagation of chemical reaction fronts [34]. Moreover, in recent years a series of papers focused on the generalizations of the original comb model, as systems encompassing superdiffusion, due to the presence of inhomogeneous convection [35] or Lévy flights [36], heterogeneous and fractional diffusion on a fractal mesh [37][38][39], and slow or ultraslow diffusion [40]. In general, comb models can be seen as a concrete representations of CTRW. A general account of these systems can be found in Ref. [41].
The analysis of the diffusion along the backbone has been obtained by performing numerical simulations of N = 3 × 10 6 RWs over a Sierpinski comb-structure, with one gasket for every backbone site. We consider RWs on several comb-structures differing in the fractal dimension of the SBs. Each SB is characterized by a primary element by means of which one can iteratively generate the fractal structure. Such a geometrical element is called δ-simplex, which is a set of δ sites joined by edges. The Euclidean dimension d of the space in which the gasket is embedded is related to the simplex by d = δ − 1 [42]. The spectral dimension d s of such structures is related to the Euclidean dimension by the relationship [42,43] while the fractal dimension d f depends on d via the relationship The behavior of the mean-square displacement (MSD) of the RW over comb-like structures is anomalous, Eq.(1), with an exponent ν depending on the spectral dimension of the SB: Of course, the unavoidable finiteness of the linear size L of the SB's makes the anomalous diffusion a transient behavior occurring before the onset of a standard regime. However, upon taking L to be sufficiently long, with respect to the mean free path, the transient anomalous behavior can be made arbitrarily long-lived. Equation (11) can be explained by a simple phenomenological argument [44] which is based on the homogenization time, t * (L), meant as the typical timescale after which the diffusion along the backbone becomes standard where D(L) is the effective diffusion coefficient depending on the scale L. For finite-size SB's indeed, the anomalous regime in the longitudinal diffusion is only transient, and soon or later, it will be replaced by the standard diffusion. The homogenization time, t * (L), can be identified with the typical time taken by the walker to "span" a single SB of linear size L. The scaling of t * (L) with L, can be easily extracted from the diffusion on fractal structures, x 2 ⊥ (t) ∼ t 2/dw , where d w indicates the random walk dimension [41]. Therefore we expect that, Once such a scaling is known, we can apply a "matching argument" to derive the exponent ν in Eq.(1), by requiring that the power-law and the linear behavior have to match at time t t * (L), We need to determine the scale dependence of the effective diffusion coefficient D(L) for the RW along the backbone. This scale dependence is d f being the fractal dimension of each SB. Indeed, when the diffusion along the backbone has reached the standard regime, it satisfies reduced by a factor f (L), which is the probability for a walker to occupy a backbone site, since only the fraction f (L) of RWs on the backbone actually contribute to the diffusion. We can safely assume that the dynamics in the homogenization regime equally visits all the sites on each SB, which become equiprobable. Accordingly, f (L) follows from a simple geometrical counting that is evident when referring to the simple onedimensional comb with finite SB of length L, [top panel of Fig.2(A)]. Clearly, each linear SB contains a number of sites N SB (L) = L/ 0 , where 0 is the lattice spacing that, without loss of generality can be set to 0 = 1. Therefore, f (L) amounts to the following simple counting: only one backbone site over N SB (L) = L total SB sites , implying that f  (14) indicates that D(L) decays by enhancing the trapping power of the SBs, this occurs for two reasons: by increasing their linear size L, or by increasing their geometrical complexity, characterized by d f . Now, using t * (L) ∼ L dw , the matching Eq.(13) yields 2νd w = −d f + d w , from which we obtain Eq.(11), providing that d w = 2d f /d s , see Ref. [16]. A similar reasoning has been generalized to compute anomalous exponents for RWs on fractal brushes, in Ref. [45].

A. Comb and comb-Sierpinski structures
Because of the coupling between the SBs and the backbone dynamics, the transport along the backbone of a comb, see the top of Fig.2(A), becomes anomalous with an exponent given by Eq. (11). Figure 3 reports the power-law behavior of the MSD along with its prediction (dashed lines) for the structures corresponding to different d s . The long-time PDFs of the RW dynamics on the backbone present the collapse onto a F ν (z) predicted by Eq.(2), upon rescaling, z = |x|σ(t), with Fig.4. To test whether the PDFs exhibit Fisher's tails, we performed a global fitting of all the curves using Eq.(5) with a ν and α as free parameters. This unconstrained or "blind" fitting procedure shows that the stretched-exponential form fits only the distribution bulk, i.e. small values of z. Moreover, the α values from the fitting deviate from the expected values, Eq.(7), see Table I. This analysis suggests that without a proper guess, there is no chance to prove that the backbone diffusion in comb-like structures follows a Fisher-like distribution, at least on the tails. An analytical prediction of the PDF can be performed using the fractional Fokker-Planck equation (FFPE) which governs  the diffusion on the comb backbone; see Ref. [46].
To start with, let us briefly recall the theory developed for the comb system in Refs. [47][48][49]. According to d the authors, the Smoluchowski equation for the particle distribution P (x, y; t) on the comb (backbone plus SB's) is The presence of the δ(y) function ensures that the diffusion along the backbone takes place only at y = 0. Two different diffusion coefficients are introduced: one along the comb teeth (D y ) and one along the backbone (D x ), with different physical dimensions. The equation for the particles diffusing along the backbone, can be derived by manipulating Eq.(15) [48] defined as where K 1/2 = D x /(2 D y ) and 0 D α t is the Riemann-Liouville fractional operator of order α [46,50]: As anticipated, Eq.(16) is a type of FFPE [46], and its solution can be compactly expressed by a particular case of the Fox function [46], also called the M-function [51,52]: Here, the MSD takes the form but most importantly, in the limit |x|/σ 1, the Fox function in Eq.(18) admits the asymptotic expansion [53,54] The expected Fisher's stretched-exponential behavior is indeed recovered, and remarkably, it reproduces the tails (z 4) of the numerical PDFs, reported in Fig.5a, without any fitting procedure. Extending this approach to the fractal combs of Fig.2 is out of the question, because the analogous of Eq.(15) cannot be drawn. However, we can pass through the CTRW representation. Let us first examine the comb model. If we consider only the dynamics along the backbone, the time spent at a site x is the return time to x, after the walker has visited the matching SB. In other words, we can imagine the particle still being at site x, while performing its Brownian motion along the finger (y). Within a certain degree of accuracy, we can assume that this residence (or waiting) time at x corresponds to the walker's first return time to the backbone (y = 0). Hence the residence time PDF ω(t) scales as ω(t) ∼ (τ /t) 3/2 , where τ is an arbitrary time-scale. On the other hand, the hopping dynamics along the backbone would ensure that the jump distribution is well approximated by a Gaussian, i.e., λ(x) = (4πσ 2 ) −1/2 exp[−x 2 /(4σ 2 )]. Following the derivation furnished in Ref. [46], the Laplace and Fourier transforms of the waiting time and jump length PDFs read respectively. Therefore the Fourier-Laplace transform of the CTRW PDF is so that the CTRW approach yields the correct MSD (19) and, most importantly, it confirms that the corresponding Fokker-Planck equation is the FFPE (15) [46]. Now let us turn to the fractal comb. The CTRW picture fully applies also to this case, with the only difference that a particle is imagined to reside at site x while wandering through the fractal hinged on the backbone's site x by its origin [see Fig.2(A)]. Hence the residence (or waiting) time corresponds to the walker's first return time to the origin of a Sierpinski gasket of spectral dimension d s , whose PDF behaves as [42] On the other side, the jump-length distribution is still a Gaussian. Thus, in this case, we have and the Fourier-Laplace transform of the walker's PDF is Inverting, the proper Fokker-Planck equation for such a process reads [46] with the MSD As for Eq. (16), the solution of Eq.(23) is again an Mfunction [51,52] satisfying the asymptotic behavior which reproduces the Fisher's stretched exponential. In general, it can be shown that the M-function appearing in Eq.(25) is related to the Lévy stable distribution with parameter ν by the following equality [55] Figure 5 shows that the asymptotic formula (26) is an excellent approximation of the distribution tails. Once again, we stress that no fit to numerical data has been implemented.
To sum up, expressions (20) and (26) clearly demonstrate that Fisher's scenario is fulfilled if we take into consideration only the tails of the PDFs (z 4), and that it is a consequence of the type of fractional equation governing the diffusion on the comb backbone.

III. RANDOM-SHEAR MODEL
Random-shear models have been proposed to study the hydrodynamic transport of solute particles in stratified porous media, under the assumption that advection is parallel to the stratification planes [22]. Let P (x, y, t) be the concentration of the solute particles and U = (U (y), 0) be the shear parallel to x depending only on the stratification height y, we can write the conventional advection-diffusion equation where the incompressibility div(U) = 0 has been taken into account and D x , D y denote the microscopic (molecular) diffusivity along x, y directions respectively. In the following, we set D y = D 0 , because we will assume that D x 0, so that the particle dispersion along x is mainly ruled by the statistical properties of U (y). In a famous paper, Matheron and de-Marsily [22] studied the anomalous transport of particles driven by an array of horizontal random-width layers with a constant velocity. The particles also undergo a Brownian diffusion along the transverse direction, therefore, they cross different layers, visiting them several times. Specifically, Matheron and de-Marsily identified two properties of the velocity autocorrelation U (y)U (0) leading to a longitudinal anomalous super-diffusion, where the average is computed over the random field realizations.
The problem (27) can be reformulated in terms of Langevin equations for an ensemble of independent particlesẋ {ξ i (t)} are independent, delta-correlated, and zero-mean Gaussian noises. Periodic boundary conditions are enforced on the y-direction at y = ±L/2, to implement a channel-like geometry along the x-axis.
Our shear longitudinal field is generated by a superpo-sition of M sinusoidal waves, i.e.
where the sum on k runs over the set k = 2π/L(1, 2, . . . , M ), with λ = 2π/L and Λ = M λ. This definition is not too restrictive for M large, since any smooth-enough field can always be expanded in Fourier modes. The amplitudes are assumed to depend on the wave vectors as U 0 is a dimensional factor that can be set to unity by a simple time redefinition, while the presence of the phases {φ k } is strongly necessary to confer a certain degree of heterogeneity (randomization) to the field. The parameter γ, defining the spectral properties of the modal decomposition of U (y), is the only quantity that will be varied in our analysis. In what follows, it will be convenient to express U (y) in the complex form For reasons of convenience that will be soon clear, we take the shear field to be anti-symmetric about the middle of the channel, y = 0, i.e. U (−y) = −U (y). This condition also sets up the phase choice to φ k = 0 or π with probability (1/2, 1/2), respectively.
The model (28,29) can be exactly solved as the equation for y is independent of x: the solution is y(t) = y 0 + √ 2D 0 w t , where w t indicates a Wiener's process, i.e. w t = 0 and w s w t = |t − s|. A substitution into the first equation yields Thanks to Eq.(30), Eq.(33) expands to Now, from expression (34), all the moments can be computed. The crucial point is to establish the meaning of the average · · · . As a matter of fact, three types of average can be carried out on this system: 1. over the noise realizations, i.e. on the Wiener's process w t : . . . w 2. over the initial conditions y 0 along the stripe: . . . 0 = L/2 −L/2 dy 0 ρ(y 0 ) 3. over the possible configurations of the disordered field U (y), i.e. over the phases: . . . φ = λ k=λ π −π dφ k /(2π). Unlike previous works on the random-shear model [22,29,41,[56][57][58][59][60][61][62][63][64][65][66], in our analysis we will not consider the average over random field realizations (phases in our case), taking only the average on the Wiener process and on the initial conditions, as done in Ref. [67], in symbols · · · w 0 . We will be assuming an initial uniform distribution of particles along the channel: ρ(y 0 ) = 1/L.
First we can compute the mean displacement (drift), which is rigorously zero because exp(iky 0 ) 0 = δ k,0 and we have used the property of w t such that exp(ikw t ) w = exp(−D 0 k 2 t).
The mean square displacement reads recalling that V k = U k e φ k /(2iM ), we have |V k | 2 = U 2 k /(4M 2 ), so that the MSD turns out to be independent of the phase disorder. The first contribution is the well-known Taylor term [68] of the effective standard diffusion, whereas the second term is the contribution leading to the anomalous transient behavior whose duration increases with the width (transversal size) L of the channel.
Before showing analytically how the anomalous regime emerges from Eq.(36), we repeat the matching argument of sec.II, in the context of the random-shear model to obtain its anomalous exponent in an intuitive way. Let us consider the case of large but finite L; of course, when t t * = L 2 /D 0 , the transverse free diffusion of particles feels the boundary effects, and the longitudinal diffusion changes from anomalous to standard. Therefore, even for the random-shear, we can assume that Eq.(12) holds true, with the exponents determined by the random shear problem. In this case, D(L) scales as according to the Taylor formula (37), where the definition, Eq.(31), has been used along with an implicit passage from the summation on k to the integral over dk.
The anomalous and the standard regimes need to match each other at time t * (L), which is the typical time, t * (L) ∼ L 2 /D 0 , taken by transverse free diffusion [Eq. (29)] to become almost uniform over the channel width L. Thus, from the matching condition, t * (L) 2ν ∼ D eff (L)t * (L), one obtains the relation (L 2 ) 2ν ∼ L 1−γ L 2 , which, by equating the exponents, leads to Therefore, as in the case of comb structures, a simple scaling argument links the anomalous exponent ν to the parameter γ defining the velocity field spectrum through Eq.(31), in addition Eq. (7) provides the value for the Fisher's parameter.
The above derivation can be made rigorous, again assuming large channel widths L 1, so the summation over k can be replaced by an integral over the interval [λ, Λ], casting Eq.(36) in the following form: where Eq.(31) has been used. Appendix A shows that, in the time interval where the constants are: Equation (41) prescribes that a particle driven by the random shear undergoes a transient super-diffusion with the exponent (38). Of course, for t L 2 /D 0 the MSD turns into the Gaussian regime characterized by a coefficient (37).
The shear parameters are kept fixed to: U 0 = 1, L = 100, and M = 100, while γ in Eq. (31) is varied in the range −1 ≤ γ ≤ 1. Moreover, our specific choice of the disorder (φ k = 0 or π) corresponds to assign random signs to the amplitudes U k = ±U 0 |k| γ/2 . The numerical integration of Eqs. (28,29) has been performed with a simple Euler scheme with a time step h = 0.01.
The first test of the correctness of the Euler integration is the agreement of the MSD obtained by the simulations with the analytical result, Eq.(36), as verified in Fig.6 by the coincidence of the red curves representing Eq.(36) with the numerical MSD (open circles). We computed also the time behavior of the moments M m (t) with m = 4, 6 to verify their expected intermediate anomalous behavior (4) occurring in the interval L 2 /(M 2 D 0 ) t L 2 /D 0 , which is represented by the dashed lines in Fig.6. For times large enough, t 2×10 4 , the Gaussian scaling of moments is recovered, with the normal exponent ν = 1/2, indicating that the system has attained the standard regime.
The insets in Fig.6 show also the alignment of moments upon raising them to the right power, [M 2m (t)] 1/m ∼ M 2 (t) (Fig.6). This alignment is a consequence of the scaling (2), and the data in the insets represent a first robust numerical test that this scaling is verified by the simulations. The same conclusion was reached in Refs. [41,56] using the average over the disordered convection fields · · · φ .
The collapse of moments in Fig.6 shows the not strongly anomalous character of the superdiffusive regime [20]. In view of Eq.(3) indeed, we have ν(m) ≡ ν = (3 − γ)/4, as is shown in Fig.7 for the same values of γ displayed in Fig.6. Now, as the scaling of the moments corroborates the validity of Eq.(2), we need to establish whether the tails of the PDFs satisfy Eq.(5) with the expected exponent α = 4/(1 + γ). Figure 8 displays the collapse of the PDF at different times according to Eq.(2). Moreover, the red dashed curves represent the fitting with Eq.(5), with the constraint α = 4/(1 + γ) being priorly imposed. In practice, only the amplitude, A, and the parameter, a, are adjusted by the fitting procedure. Although the fitting curve fails to reproduce the bulk of the simulated PDF, it is very reasonable for the far tails. Besides the values γ = 0, 0.4, the consistency of the Fisher's scenario has been also verified for γ = 0.2, 0.6, 0.8 (not shown).
However, the Fisher's prediction seems to fail for the PDF with negative γ; see the red dashed curve in the bottom panel of Fig.8.Although the tails of the rescaled distributions are still stretched exponential, the exponent is different from α = 4/(1 + γ).
As a final remark, we can say that the behavior of the longitudinal PDF is strongly dependent on the properties of the shear field, so the Fisher scenario is not so robust. For instance, if the phases are such that φ = ±π/2, with probability 1/2, we obtain a distribution that still follows the scaling law (2). However the tails are not Fisher-like (see Fig.9) also because a generic disorder does not preserve the symmetry, x → −x, of the Fisher's distributions.

IV. CONCLUSIONS
In this work we have studied the subdiffusion along the backbone of random walks on comb-like fractal structures and the superdiffusion of the Lagrangian dynamics of Brownian particles in a random-shear velocity field. In this case the anomalous transport occurs along the shear direction.
In the comb systems, the anomalous exponent can be analytically obtained as a function of the fractal dimension of the Sierpinski sidebranches decorating the backbone. In the shear model the anomalous exponent is a function of the parameter γ defining the spectral properties of the Fourier-mode combination of the velocity field.
We focused our analysis mainly on the scaling properties of the spatial PDF of the process along the preferential transport direction: the backbone (for combs) and x-axis (for random shear), because we were interested in testing the assumption that such PDFs develop Fisher's stretched-exponential tails exp(−a|z| α ), with α related to the anomalous exponent by the relation α = 1/(1 − ν).
Our simulations show that for comb systems, the PDF of the anomalous sub-diffusion along the backbone follows the Fisher's tails, regardless of the complexity of the fractal sidebranches. This numerical finding has been also supported by analytical predictions based on the extension of the fractional Fokker-Planck equation to comb fractal structures.
For the random shear, we found that the transport is anomalous with an exponent ν = (3 − γ)/4. The particle PDFs along the x-axis exhibit stretched-exponential tails with an exponent consistent with the Fisher's prediction, Eq.(39), only for γ > 0. Surprisingly, for γ < 0, numerical simulations seem to indicate that the Fisher's scenario breaks down.  γ). Notice that the case γ = 0, for which ν = 3/4, corresponds to the Matheron de-Marsily model [22]. The last panel, γ = −0.4, suggests that the Fisher's fitting is not satisfactory for negative γ.  2), is verified, however, the PDF, both for symmetry and tails, is not in the Fisher's class.