Transport in a stochastic double diffusivity model

A recent study analyzed the role of stochastic perturbations on the interface dynamics of two interacting species within a double-diffusivity framework, involving double diffusion models. The model relied on a restricted translation–reflection (TR) symmetry manifold, leading to a single variable description. The present study generalizes this model for a TR symmetry violating system that does not permit reduction to a single variable dynamics, leading to a hitherto unseen stochastic resonance (SR), a mechanism that indicates discrete, rather than a continuous, mode of energy transport. The SR exhibited by the model captures the signature fast transport observed in stochastically driven dynamics of nanopolycrystals, that previous deterministic models failed to emulate. We show that the speed of transfer relates to the strength of energy cross-correlations between the two diffusing species communicating through interface dynamics that eventually drives the energy throughput and identifies the role of stochasticity in nanopolycrystalline transport.


Introduction
A continuum theory of diffusion for media containing two types of "low" and "high" diffusivity paths was initially proposed by Aifantis [1][2][3] to supplement existing discrete models for transport along dislocation lines (dislocation pipe diffusion) and grain boundaries (grain boundary diffusion). It is described by two second order "Fick type" partial differential equations for modeling the transport of the diffusive species along each of the existing paths, containing also a "chemical reaction type" or source term modeling the spontaneous jumps of diffusion species from one path to another. Uncoupling of these double diffusivity (DD) model equations resulted in a higher fourth order partial differential equation which contained new terms characterizing previously proposed transport equations, namely: a second time derivative term found in the telegrapher's equation [4]; a mixed time-second space derivative reminiscent of Barenblatt's seepage equation in fissured rocks [5]; and a fourth spatial derivative term appearing in the Cahn-Hilliard theory of spinodal decomposition [6]. The DD model was also employed by Aris (in a landmark paper dedicated to the memory of Dick Bellman) who also extended it to discuss diffusion in membranes with n pathways. Details pertaining to the original derivation of the DD model can be found in [7], where a "mechanical basis" for interstitial diffusion was provided in terms of the momentum balance for the diffusing species and the introduction of a "drag" or "diffusive" force to account for the momentum exchange between the solute and the solid matrix (in analogy to Maxwell's statistical mechanics theory of polyatomic gases). Mathematical properties of the DD model and analytical solutions expressed in terms of known fundamental solutions of the classical diffusion equation were provided in [8][9][10].
The DD model has found numerous applications in a variety of problems and configurations where two distinctly different paths for transport (lattice and dislocation pipe or grain boundary diffusion; heat or electricity conduction in inhomogeneous media; fluid or traffic flow; and others) are present [11][12][13][14][15]. Within a continuum formulation, this is equivalent to assuming that two non-equilibrated field variables (two concentrations, two temperatures, two pore pressures) co-exist in the elementary volume at hand: e.g., "bulk" vs. "grain boundary" solute concentration, "electron" vs. "ion" temperature; "pore" vs. "fissure" pressure. Our present work focuses on diffusion in nanopolycrystals (where the "grain boundary" space occupies 30% of the "bulk"), within the context of a stochastic extension of the DD model. More specifically, we elaborate on the coupled system of the two partial differential solved in [9,10] and employed in [16] to interpret experimental measurements of diffusion in polycrystals and nanocrystalline aggregates.
As reviewed in [1][2][3], the first models for mass transport in high-diffusivity paths (e.g., dislocations or grain boundaries) were discrete in nature, assuming a single dislocation with a core size R 0 or a grain boundary with thickness H 0 placed in a medium of infinite extent and solving the classical Fick's equation for these configurations, with appropriate boundary conditions, to determine the concentration-distance profile along the dislocation "pipe" or the grain boundary layer. In contrast, the continuum DD model assumes a continuous distribution of high-diffusivity paths (with volume fraction "f" as in "mixture theories" for composites) and allows for both mass transport and mass exchange between the two distinct families of locally co-existing paths within the elementary material volume considered. The continuum DD model becomes increasingly relevant with decreasing grain size, where deterministic (and more importantly stochastic) heterogeneity effects (dislocation clusters, grain boundary ledges, etc.) prevail over the idealized geometry used in the discrete models. In fact, with decreasing grain size down to the nanometer level, dislocation core sizes and grain boundary thicknesses are safer to be represented by the aforementioned average local volume fraction parameter. This, in turn, leads to the definition of two average local concentration fields, each obeying Fick's law with distinct diffusivities for each family of paths and also accounting for mass transfer between them. The dependence of granular flow on rheology and consequent impact of local friction has also been shown to abide structural invariance in plasticity models [17]. Apart from the obvious importance of such observations from the perspective of phase transitions and implicit renormalization of the diffusion constant, presence of multiscale features as a consequence of microstructure evolution is a well-studied phenomenon in metal plasticity, particularly in the understanding of the impact of grain boundaries [18]. Topical experiments on Copper, using optical and electron microscopy along with micro x-ray tomography, to analyze the nature and degree of damaged microstructure show tensile damage in the recovered samples is examined using [19]. This latter study points to differential impact of grain sizes on surface velocity, another feature that has been traditionally beyond the scope of traditional DD models, an excellent summary of which has been outlined in [20]. As this article [20] explains, dislocation patterning in thin films was already understood to be a multiscale feature that could be analyzed by gradientdependent modeling, but micro-nano level stochastic fluctuations remained outside the scope of such studies. A key study analyzing the nature of distribution of such granular imperfections and the consequences of that on the emergent statistics of such processes was studied recently [21]. The present study adds the critical element of quantifying the dynamical evolution of such nanolevel granular processes by analyzing domain boundary movement using Langevin and Fokker-Planck based models.
Moreover, for the important class of ultrafine grain (UFG) materials (grain size of ∼100 nm fabricated by severe plastic deformation/SPD), a new family of high-angle (non-equilibrium) grain boundaries emerge with different substructure than their usual low-angle (equilibrium) counterparts and, thus, different diffusion properties. The same is true for other types of novel nanocrystalline (NC) materials (grain size of ∼ 10 nm fabricated by vapor deposition) where the density of volume fraction of triple grain boundary junctions (TJs) increases significantly, leading to another distinct family of paths with very high diffusivity, thus necessitating the use of a "triple" diffusion model [22]. The multiscale and nonlocal nature of diffusion at the nanoscale leads to hybrid or anomalous transport of deformation energy, characterized by the coexistence of a diffusive field with a propagative component [23]. Such anomalous transport has been shown recently to be captured by fractional diffusion equations, able of modeling both retarding subdiffusion as well as accelerating superdiffusion [24]. Such anomalous diffusion possessing non-unique diffusion coefficients (diffusion coefficients varying with time), corresponding to cases with rapidly-changing microstructure are not within the scope of the present formulation which only considers two distinct diffusion paths as well as small random perturbations of the materials' microstructure modeled through the noise terms entering the diffusion equations.
More detalis on the physics and properties of high density or volume fraction grain boundary materials can be found in the materials science literature (e.g., [25,26]). For micropolycrystals (grain size ∼ 1 μm), the two families of paths correspond to the bulk (grain interior) and the surrounding internal surface/ grain boundary space. In all these configurations, the applicability of the DD model renders new insight for diffusion in inhomogeneous media and provides new possibilities for interpreting related experimental results, as is further demonstrated in the present article.
The above discussion clearly points to the lack of appropriate consideration of internally induced fluctuations arising out of the structural difference between the two types of diffusion paths involved, exhibiting different local material topology and structural randomness. In line with well established stochastically forced flow models [27][28][29][30] representing archaic dynamical randomness (generated close to the boundary layers), the present article will explore this realistic limit of double-diffusion, thereby accounting for all modes of randomness. It will be structured within a well-knit Langevin formulation of stochastic dynamics [31,32]. Phenomenologically, this can be viewed as the effect of an external stochastic force that is randomly redistributing the relevant spatial material inhomogeneity: for example, that of the high-diffusivity paths, in which a stochastic increase (decrease) in temperature or fluctuation of internal stress, may extend (or contract) the interlayer grain boundary thickness between two nanosized grains and, thus, alter the configuration of the associated structural defects (dislocations, disclinations). Such a multi ensembled stochastic reorganization of the local space distribution is sometimes known to create new universality classes [32,33], another aspect that could be of critical interest in analyzing the dynamics of nanopolycrystals.

The model equations
The double diffusivity mathematical models referred to in the previous section were purely deterministic and the role of stochastic fluctuations, which become increasingly important as the specimen size decreases, was not considered. The first relevant work on the DD model under stochastic forcing, analyzing the implications of micro/nano structure-induced concentration fluctuations, was recently done by Chattopadhyay and Aifantis [34]. In a major departure from the general trend of results derived from deterministic models, this work highlighted the role played by stochastic effects in stretching the relaxation time scales of nanopolycrystal samples. Moreover, the implications of boundary layer perturbations on inhomogeneity induced micro-and nano-fluctuations were analyzed. What this study did not reveal, though, was the quantitative importance of the nature of noise distribution that drives the force field. Also, the role of multiple time scales could not be unambiguously dealt with, as often observed in experiments with such nano-sized samples.
As established in [34], we focus here on the time evolution of the concentration of two stochastically interacting diffusing species, identified in the "slow" ( 1 ) and "high" ( 2 ) diffusivity paths. Assuming a minimalist view for the interaction between these two diffusing populations along two different paths (e.g., intracrystalline grain interior (GI) and grain boundary (GB) or intercrystalline (IC) and triple grain boundary junction (TJ) space) with corresponding concentrations 1 and 2 , we can represent the dynamics by an extended stochastically forced DD model as follows: are Ornstein-Uhlenbeck noises [35] that are introduced to account for substructural randomness associated with both boundary layers and the surrounding bulk matrix. They are defined, as usual, by the relations where � i (0) denotes noise strength and is the relaxation time, is the usual Dirac delta function, and ij is the Kronecker symbol (no summation of the dummy indices implied). In Eqs. (1a) and (1b) above, D 1 and D 2 respectively refer to the diffusion coefficients in the two paths, with the constants 1 and 2 denoting the respective formation or depletion rates of concentrations of the diffusing species which can pass from one path to another. Equations (1a) and (1b) represent interacting diffusion fluxes, incorporating linear interaction in between the two diffusive paths. For all practical purposes, we will resort to the non-dimensional representation of this model, defined as follows: Here is (x, t) (i = 1,2) are the scaled variables derived from 1 and 2 , respectively, such that , c 3 = 1 , d 1 = d 2 = 1 and d 3 = 1∕c 2 , whereas 1 and 2 are the scaled stochastic forcing terms, such that (1a) Here i (0) ( i = 1, 2 ) are the rescaled noise strengths corresponding to 1 and 2 respectively, and is the noise decay time. Our primary focus in this work will be to analyze the role of the Ornstein-Uhlenbeck decay time in determining the scaled diffusing species in the dynamics of the IC-TJ paths.

Phase evolution dynamics
As detailed in [34], the phase evolution kinetics of the IC and TJ diffusion will be discussed using the standard language of two-point correlation functions involving the scaled variables. Henceforth, without loss of generality, we will identify the variables in Eqs. (3a) and (3b) as 1 (x, t) and 2 (x, t) . Fourier transformation of these variables, can be represented where generically represents either i or i (i = 1, 2). This will lead to the following set of coupled linear equations where the matrix M is defined as To analyze the kinetics from this stochastic model, we need to calculate measurable quantities. This is done using the formulation discussed in [28,30]; the quantities that could be experimentally measured will be the (Brownian) root-mean-square averaged quantities of their deterministic equivalents as detailed below: and a complementary set of cross-correlation functions given by 12s The above cross-correlation term shown in Eq. (7) is likely to be much smaller compared to the root-meansquare (rms) quantities defined in Eq. (6) above. The quantities within the brackets ("< . > ") indicate "ensemble average" while * is is the complex conjugate of is .

Phase autocorrelation and cross-correlation
The kinetics of our stochastic model is estimated using the language of two-point structure functions, popularly referred to as autocorrelation and spatiotemporal functions, separately between similar variables as well as between different ones (crosscorrelation). Section 4.1 below shows the autocorrelation structure, Sect. 4.2 refers to the cross-correlation between zero-point different autocorrelators. The following Sects. 4.3 and 4.4, respectively, allude to spatial and temporal correlation functions, both for like minded variables as also for cross-correlations.

Autocorrelation
To evaluate the scaled is rms ( i = 1, 2 ), we need to calculate their corresponding complex conjugates in the Fourier transformed k − space represented as follows It is important to note that unlike the case with Gaussian white noises [34], the presence of Ornstein-Uhlenbeck dynamics, leads to a non-trivial spatiotemporal correlation for the noise part due to finite time correlation of fluctuations. This can be estimated to give the following form VT , V and T being the total volume and time of the evolution. In the k − space, the complex correlation functions can be calculated to obtain the following forms . The denominator 1 is given by . The equation 1 = 0 defines the poles of the scaled model. The removable singularities lead to the following 4th-order polynomial pole structure ( = Ω are the poles): The correlation functions shown in Eqs. (10a, 10b) can now be used to evaluate (analytical results validated against Mathematica 12) with their corresponding plots given in Fig. 1, from where it can be seen that the autocorrelation function auto 1s shows much higher saturation rates, but is three orders of magnitude lower than auto 2s . (13b) As can be easily seen from Eq. (13a), for Ω = 1 , the energy associated with either energy will diverge. In experimental parlance, this is the point of resonance, mediated by the presence of stochastic noise, and hence can be identified as "stochastic resonance" (SR) [36]. As shown in a complementary model [37], this is a unique result as conventional SR is known to be a nonlinear effect only, whereas here the dynamics is entirely linear. This resonance point also gives us a quantitative idea of the validity regime of the Ornstein-Uhlenbeck noise strength for this dynamics that comes out as ≤ 1 . This -inequality sets out a validity regime for the lower limit of the wave-vector k = k 0 given by not a particularly interesting case for nanopolycrystal double diffusion, as well as for situations involving diffusion of two different species through the same media (obtained by carefully choosing values of the other parameters), this can serve as an appropriate benchmark against which relative growth of the 1 or 2 species could be measured. Comparing with Eqs. (12a, 12b), the stochastic resonance conditions implies that while very high diffusivity rates D 1 and D 2 are not allowed for such a resonance (as Ω s will be imaginary), large values of "mixing rates" k 1 and k 2 ensuring → 0 will precipitate stochastic resonance. In other words, the impact of heterogeneity will be expressed through the mixing rates and not diffusivity itself due to the specific temporal scale chosen by this type of resonance. This feature is unique compared to conventional "internal length" theories. that Here the average energy dissipation rates of each of the two phases (IC and TJ equivalent) will be equal to each other. For many experimental observations (e.g., [16]) concerning double diffusivity, the two diffusive constants typically differ by about three orders of magnitude whose exact correlation forms can be analyzed using Eqs. (13a,13b). Another key qualitative difference between the two cases D 1 = D 2 and D 1 ≠ D 2 is the lack of "stochastic resonance" in the former.

Cross-correlation
While diffusive relaxation of individual species can be measured against the autocorrelation thresholds 1s auto or 2s auto , the off-diagonal terms in the Hessian will contribute to a diffusive flux across different species that, especially when long-ranged, could contribute substantially to the transport of either species. The scaled cross-correlated autocorrelation function 12s is given as follows (analytical results validated against Mathematica 12) As Fig. 2 clearly indicates, the qualitative nature of transport remains the same as that of individual diffusive relaxation as depicted in Fig. 1, with the flux magnitude for cross-correlated diffusion being of the (14) [ 12s (auto-cross)] 2 same order of magnitude as auto 2s . The trend remains unchanged for = 0.03 (dot-dashed), 0.04 (dashed) and 0.05 (dotted), although there are minor quantitative differences in their saturation rates. Since Figs. 1b and 2 are conjugate of each other for identical noise strengths, e.g., 1 = 2 = 1 , as represented by equations (13b) and (14) respectively, we chose complementary pairs for the respective figures; in particular 1 = 5, 2 = 0.1 for Fig. 1b and 1 = 0.1, 2 = 5 for Fig. 2.

Spatial correlation of phases
In this section, we analyze how the concentration of each phase changes with spatial distance in the dynamical equilibrium limit. This can be done by solving the respective spatial correlation functions of each phase for all times and then taking ensemble averages over all noise realizations, an approach in line with [42]. By definition, we have This then leads to the following (radial) spatial correlation functions Within a very short interval, contributions to the correlation functions from such large-spatial separations can be seen to decay to a saturated state as shown in Fig. 3. This confirms the convergence of the integrals in Eq. (16a). In Fig. 3, we show how the spatial correlation oscillates and then saturates with increasing separation distance, for the case 1 = 2 . The oscillation relates to the phase contribution that numerically corresponds to the imaginary part of the solution; the decaying envelope confirms that spatial correlation trails off as the separation between the active sites increase.
As with autocorrelation, the cross-correlation of the two diffusivity paths will lead to a source-sink (or reverse) effect between the two paths 1s and 2s , (16b) as depicted in Fig. 4. The corresponding correlation function is given below: There is a remarkable quantitative similarity between cross-relations in the autocorrelation and spatial correlation sectors. As Fig. 4 shows, the magnitude of spatial cross-correlation is almost comparable to that of one of the individual spatial two-point correlations and three orders of magnitude larger than the other.

Temporal correlation of phases
It is well known that a fundamental consideration in multi-phase systems is the time evolution of the interface separating two different phases. Often such systems are known to be stochastically perturbed and hence non-equilibrated in nature, potentially (17) [ 12s (rms-cross, rendering the relevant dynamics as oscillatory, or with oscillatory-rotatory instability leading to chaos [38]. Phase control through synchronization driving such systems away from the chaotic bifurcation point has, in fact, benchmarked the hare-lynx model in ecology [39]. By analogy with two-phase systems, we calculate below the theoretical quantities which, in principle, can be compared with the experimental set ups describing the temporal correlation dynamics. In the present case, the relevant variables necessary to model such a stochastically driven two-phase system are i rms (T) , which are encapsulated in the following equations for i = 1, 2 . The scaled temporal correlation functions are given by (18) [ is As is quite obvious, both temporal correlation functions for 1 and 2 show stochastic resonance for Ω = 1∕ . With longer time gaps between the recorded fluctuations, there is an evident decay in the strength of correlation. Figure 5a and b, show a steady decline with increasing T where the faster diffusive energy relaxation leads to stronger connection with fluctuations separated by larger time gaps. As like the auto and spatial correlations, cross-correlation between the paths would also impact the temporal correlation, defined in Eq. (20) and plotted in Fig. 6.
A look at Fig. 6 indicates that the amplitude of temporal cross-correlation is between that of the slowly and that of the quickly diffusing path, with the decay profile similar to both. The results are similar to that of the spatial cross-correlated paths. The crossover at around T = 0.1 is indicative of the time point at which the conservative deterministic dynamics is taken over by the additional energy pumping through the cross-correlation terms. As expected, a larger decay time, amounting to a longer relaxation span, is precipitated by a slower release of energy.

Discussion
In this article, we have extended the continuum double diffusivity (DD) model to include stochastic forcing, in order to analyze the transport dynamics of nanopolycrystals. While earlier works have established the origin of stochasticity [34,37], a simple hand weaving recapitulation can help in reminding us about the importance of stochastic contribution to the DD modeling architecture. The stochasticity in the model in Eq. (3) relates to a combination of boundary layer fluctuations and separate contributions due to bulk inhomogeneity in the nanopolycrystaline sample, leading to randomness in the distribution of the diffusion profiles. A key aspect of stochastic models is the "avalanche" effect whereby minor random perturbations accumulate and pile on [40,41].
Moreover, previous higher order gradient material mechanics models (including DD ones), deterministic or stochastic, initially structured on multiscale description [43], commonly rely on the existence of translation-reflection (TR) symmetry, whereas experimental description often hints at the violation of such symmetries in finite sized systems [16,25,26]. The present study sacrificed the aforementioned minimalist description in favor of a rigorous one. The violation of the TR symmetry not only leads to inhomogeneity in the definition of the underlying model, but also the transport process sees a major energy transduction through stochastic resonance. The modeling results highlighted here, using input parameters from [16], leads to key quantitative understanding of the interface dynamics of dual diffusion in nanopolycrystals. In particular, the role of heterogeneities and lattice imperfections, largely overlooked in [16], can now be analyzed using the syntax of stochasticity, again emphasizing the importance of correctly interpreting randomness in studying experimental data. Our results show that root-mean-squared auto and cross-correlations in a stochastically forced kinetic model, e.g., that for phase densities in a forced Fick's type model, actually play the role of their equivalent quantities from deterministic models, e.g., a conventional Fick's model [44]. Technically, we introduce stochasticity as a white noise in a Fick-type model [34,37,44] of coupled equations describing the concentrations at each type of diffusion paths. Comparisons between the two interactive diffusive paths (or corresponding concentrations) have been made in terms of their spatiotemporal crosscorrelations and autocorrelations. As we see, the energy transfer is aided by the cross-correlations of the two diffusivity paths, in addition to translational energy fluxes. As discussed in detail in [31], microscopic fluctuations either in the nature of forcing or as part of the dynamical evolution process of a system, typically modeled through a combination of Langevin and Fokker-Planck equations, serve as additional non-equilibrium energy sources that are not present in deterministic models. This may lead both to additional energy flux, expressed as stochastic resonance [37] in the present DD model, or as shearing stress in similar fluid flow problems [28][29][30]. In passing, we point out that local internal stress or substructure fluctuations are not accounted for in recent generalized formulations of elasticity, plasticity, transport and coupled thermo-chemo-mechanical models. As such models are now increasingly extended from macro/meso to micro/nano scales, it is advisable to enhance them with combined deterministic-stochastic terms to interpret material behavior at these scales (see, for example, [43]).
The autocorrelation plots in Fig. 1a and b conform qualitatively to predictions based on the deterministic double diffusivity models (e.g., [16]); although the amplitudes are higher due to additional energy inputs through stochastic forcing. As to the cross-correlation terms in the spatiotemporal dynamics (Figs. 3a,  b, 4), the stochastic model leads to a new regime of description where stochasticity mediates off-diagonal terms, often asymmetric forcing across multiple variables even at the first Gaussian approximation order. While this is very much an expected part of real life nanodynamic processes, the deterministic DD model does not capture this aspect that we have successfully addressed now.

Conclusions
Models of multidiffusion are nothing new, neither is the knowledge about their ability explain the dynamics of physical processes with multiple time and length scales [1-3, 9, 10]. In a series of paradigmatic works, Aifantis, et al showed that the "Internal Length" theory could unilaterally summarize such multicomponent transport. These studies initially lacked the stochastic component [8], but was integrated with later models [11,12,14,45]. What was lacking still was the ability to associate experimental observations of transport processes, specifically multispecies diffusion rates [22,46], with corresponding models. That knowledge gap got recently bridged by showing that thermally driven stochastic forces arising from inherent inhomoegenieties and lattice defects could actually contribute a lot more than previously envisaged [34,37].
A key success of these stochastically driven models [34,37] was in accurately quantifying established experimental results of multidiffusion models [16]. These studies [34,37] had a technical simplification embedded though, that of combining the two coupled Fick equations of motion to obtain a single higher ordered PDE-equation that was analyzed to explain experimental results not only for double-diffusion models, but also in transdomain studies related to calculation of infection kinetics [47]. The issue was threefold. First, higher ordered models have structural singularities that may prompt numerical artefact in predictions. Second was a more fundamental issue concerning translational and rotational symmetry violation that such a transformation (from two coupled diffusion models to a single PDE model) could amount to. Finally, and most importantly, these stochastic models although correct in their engineering remits were still unable to analyze the nature of energy transfer during transport. They could not explain that such transports are essentially discrete processes where maximum energy transfer has to happen during stochastic resonance, and not as a continuous energy transfer process. This article completes this theoretical loop of understanding.
An important aspect of this analysis is the relative independence of both spatial and temporal correlation functions to the stochastic fluctuations. Over a wide range of noise strengths ( 10 −9 < 0 i < 10 −2 , i = 1, 2), the correlation functions showed no remarkable qualitative change and only minor quantitative changes, thereby confirming the stability of this model to noise perturbations. This, indirectly, explains why some earlier theories [16] have managed to model experimental results reasonably well for some cases, while faltering in others. For both spatial and temporal cross-correlations, the magnitudes are found to be comparable with that of their slow counterparts, a direct implication of long-ranged correlation amongst "fast" and "slow" diffusion paths. The results are very likely to see drastic change if nonlinear convective effects are introduced, as in fluid mixtures. The bulk energy transfer through stochastic resonance is a unique feature of prevailing randomness, and this should be explored further on the basis of existing or future experiments.
In summary, we point out some interesting features that stochasticity brings into the double diffusivity model by forcing a conservative energy flux (i.e., rate change of the total diffusivity is zero) with stochastic perturbations. This is only a first step toward integrating the double diffusivity properties with real life fluctuations that could modulate the process. As this model and the corresponding higher-order diffusion equation have been shown to effectively interpret transport in heterogeneous media possessing more than one family of transport paths, as well as nonlinear phenomena at the nanoscale, more detailed analysis focusing on the precise nature of randomness, as also delving deeper into the microscopic dynamics of the process, need to be conducted in the future.