Continuous Time Random Walk in a velocity field: Role of domain growth, Galilei-invariant advection-diffusion, and kinetics of particle mixing

We consider the dynamics of a separable Continuous Time Random Walk (CTRW) when the random walker is biased by a velocity field in a uniformly growing domain. Concrete examples for such domains include growing biological cells or lipid vesicles, biofilms and tissues, but also macroscopic systems such as expanding aquifers during rainy periods, or the expanding universe. The considered CTRW can be subdiffusive, normal diffusive or superdiffusive, including the particular case of a L\'evy flight. We first consider the case of zero velocity field. In the subdiffusive case, we reveal an interesting time dependence of the kurtosis of the particle probability density function. In particular, for a suitable parameter choice, we find that the propagator, which is fat tailed at short times, may cross over to a Gaussian-like propagator. We subsequently incorporate the effect of the velocity field and derive a bi-fractional diffusion-advection equation encoding the time evolution of the particle distribution. We apply this equation to study the mixing kinetics of two diffusing pulses, whose peaks move towards each other under the action of velocity fields acting in opposite directions. This deterministic motion of the peaks, together with the diffusive spreading of each pulse, tends to increase particle mixing, thereby counteracting the peak separation induced by the domain growth. As a result of this competition, different regimes of mixing arise. In the case of L\'evy flights, apart from the non-mixing regime, one has two different mixing regimes in the long-time limit, depending on the exact parameter choice: In one of these regimes, mixing is mainly driven by diffusive spreading, while in the other mixing is controlled by the velocity fields acting on each pulse. Possible implications for encounter-controlled reactions in real systems are discussed.

In the original version of the CTRW model, the probability density functions (PDFs) of the waiting times (also called trapping or sojourn times) between successive jumps and of the lengths of individual jumps are assumed to decouple, that is, are independent of each other. In the case when the PDF of the waiting times τ is fat-tailed and scale-free, ψ(τ ) ≃ τ −1−α with 0 < α < 1, and the variance of the jump length PDF is finite, the long-time limit of the model is known to yield anomalous diffusion in the subdiffusive range with the mean-squared displacement x 2 (t) ≃ t α [1][2][3][22][23][24][25][26].
In the more general case, given specific forms of the waiting time and jump length PDFs, the emerging behaviour may be subdiffusive, diffusive, or superdiffusive. This versatility of the model, together with the fact that it can be shown to be equivalent to a generalised master equation [40] and a fractional diffusion equation (FDE) in the anomalous diffusion case [41][42][43] makes the CTRW a popular choice to model anomalous transport. One of the advantages of the FDE formulation is that it can incorporate the effect of external force fields and various boundary conditions in a natural and transparent way [23,24,44]. Further generalisations are also possible, for instance, accounting for the effect of finite lifetimes of tracer particles (so called "evanescent" or "mortal walkers") [7,[45][46][47], or of chemical reactions occurring at random times and locations [48][49][50][51][52][53][54][55][56][57][58].
Recently, it has been discussed how separable CTRWs need to be formulated when the domain, on which the process is running off, is itself explicitly evolving in time. In particular, it has been shown that an FDE can be derived [59][60][61][62], whence the case of normal diffusion on an evolving domain [63] is recovered in the appropriate limit. The obtained FDE applies when the diffusing particles stick to the evolving domain, implying that they experience a drift even when they do not jump. The interplay between diffusive transport and the drift associated with the growth or contraction of the embedding medium gives rise to the onset of striking effects. These include an enhanced memory of the initial condition [59,63] and the slowing-down and even the premature halt of encounter-controlled reactions [64][65][66]. In the case of subdiffusive particles evolving on an exponentially shrinking domain, a so-called Big Crunch may happen. This phenomenon was first discussed in reference [59]; it consists in the collapse of an initial particle distribution with finite extent to a delta function as a result of the strong localisation caused by the domain contraction.
Concrete examples of expanding (or shrinking) domains include biological cells in interphase [67], growing biofilms [68,69], growing biological tissues [70,71], and growing or shrinking lipid vesicles [72,73]. The latter can be controlled easily, for instance, by adjusting osmotic pressure in solution [74]. Water drops, puddles, or aqueous solutions in a Petri dish shrink simply by evaporation [75]. On a geophysical scale, groundwater aquifers may be recharged by major flood events and thus the volume for tracer dispersion increased [76,77]. On Earth, the subducted oceanic lithosphere is stretched by the convective mantle, and both are homogenized, among other mechanisms, by diffusion [78,79]. Finally, expanding domains are traditionally considered in cosmological models describing the diffusion of cosmic rays in the expanding universe [80,81].
In the present paper, we consider the case where both a domain growth process and a left-right bias of the random walk are simultaneously at play. In reference [61], such a combination was considered for the case where the bias stems from a force field that only manifests itself at the time of each jump. To model the effect of the force field, a non-symmetric jump length distribution was used, and the corresponding Fokker-Planck equation (FPE) was obtained. This approach has been recently used to deal with the Ornstein-Uhlenbeck process on a growing domain [82].
In contrast, here we will focus on the case of a bias arising from a velocity field. This field is still at play while the walker rests between jumps. Direct realisations of such a situation could occur in biological cells in the presence of active, motor-driven motion [83], in suspended giant vesicles simply by gravity. In subsurface aquifers the flow field corresponds to groundwater streams towards a spring or well.
In a static domain, a constant force field and a constant velocity field yield the same type of advection-diffusion equation as long as the particles are normal-diffusive. However, this is no longer true when the CTRW becomes subdiffusive [84]. Then, a constant force field is assumed to act only on particles when they are not trapped, while in a constant velocity field the particle is constantly advected. The former case may, for instance, correspond to charge carriers in amorphous semiconductors (in which they are trapped at impurities) in the presence of an electrical field [2] while the latter may correspond to a particle moving in a flowing complex environment such as an actin gel. Of course, this lack of equivalence carries over to the case of an evolving domain. In [85] a fractional diffusion-advection equation (FDAE) was derived for a subdiffusive CTRW in the presence of a constant velocity field. One of our main goals will be to extend their result by considering a CTRW which takes place in a uniformly growing domain, and also including the case of superdiffusion.
In our derivation of the sought FDAE, we will first study the case in absence of the velocity field, to see that the domain growth itself may induce interesting behaviour of the moments of the particle position. For instance, in the subdiffusive case, when the physical domain grows with a power-law rate we see that a fat-tailed propagator may evolve into a Gaussian-like propagator for a suitable parameter choice.
After deriving the FDAE we will use it to study the mixing kinetics of a pair of diffusive pulses evolving on a growing domain in the presence of the velocity field, and we will discuss possible implications of the results for the kinetics of encounter-controlled reactions in real systems [64][65][66].
The remainder of this paper is organised as follows. In section 2 we first recall the main results for a symmetric CTRW on a static domain, and we subsequently discuss how the kurtosis of a subdiffusive walk changes when the initial domain grows uniformly in time. In section 3 we carry out a similar programme for a CTRW subject to the action of a velocity field, and we will derive the relevant FPE for the case of a uniformly growing domain. In section 4 we study the mixing of diffusive pulses that are biased by velocity fields acting in opposite directions. Finally, in section 5 we summarise our main results, discuss their possible relevance for encounter-controlled reactions in real systems, and outline possible extensions of the present work.

CTRW in the absence of the velocity field: Static versus growing domain
In this section we compare the behaviour of a separable, symmetric CTRW evolving on a one-dimensional static domain with the same walk on a uniformly growing domain. For both cases we discuss the difference in the behaviour of the moments of the particle position. In the case of the growing domain special emphasis is paid on the interesting onset of a time dependence at the level of the kurtosis. This section is also intended to introduce the general concepts and thus prepare the reader for the case of a CTRW subject to a constant velocity field, discussed in section 3.

Static domain
We start with a brief reminder of the derivation of a bifractional equation describing the diffusive limit of a fat-tailed CTRW. For further details we refer to the review [23] and to the recent monograph [48].
Consider a particle performing a one-dimensional, symmetric CTRW with decoupled jump length and waiting time PDFs, respectively denoted by λ(y) and ϕ(t).
Since the random walk is symmetric, the jump length PDF reflects this symmetry, λ(y) = λ(−y). The Fourier-Laplace transform W 0 (k, u) of the particle's position PDF W 0 (x, t) is known to obey the Montroll-Weiss relation [1][2][3] where W 0 (k, 0) = W 0 (k, t = 0) denotes the (Fourier-transformed) initial condition, ϕ(u) is the Laplace transform of the waiting time PDF, and Φ(u) = u −1 (1 − ϕ(u)) is the Laplace transform of the sticking probability Φ(t) = 1 − t 0 dt ′ ϕ(t ′ ) for not performing a jump up to time t. Here, we have used standard definitions of the Laplace transform and of the Fourier transform The subscript "0" in the definition of W 0 indicates that this quantity refers to the case where the velocity field is absent. The diffusive limit of the CTRW process and the associated FDE are obtained from the long-time behaviour of ϕ(t) and from the large-|y| behaviour of λ(y), which respectively correspond to the small-u behaviour of ϕ(u) and to the small-k behaviour of λ(k). For these transforms of the respective PDFs we use the well-known forms [29] with 0 < α ≤ 1, and with 0 < µ ≤ 2. This implies the asymptotic power-law forms ϕ(t) ∼ τ α /t 1+α for t ≫ τ with 0 < α < 1 as well as λ(y) ∼ σ −µ |y| −1−µ for |y| ≫ σ with 0 < µ < 2.
The GL fractional operator has the property [86] This fractional operator is equivalent to the Riemann-Liouville (RL) fractional derivative provided that both operators are applied to sufficiently smooth functions f (t) at t = 0, that is, when the condition [86] holds. Unless otherwise specified, we will henceforth assume that this is the case and therefore use the RL fractional derivative in what follows. The RL operator 0 D 1−α t is simply the first derivative of the RL fractional integral For the special case µ = 2 implying a finite variance of the jump-length PDF, it is possible to obtain the behaviour of the moments associated with the solution W 0 (y, t) of equation (8) either from the exact solution or from the corresponding hierarchy of differential equations. For 0 < µ < 2, in contrast, only fractional moments of order ν exist as long as 0 < ν < µ. For a particle initially located at the origin, W 0 (y, 0) = δ(y), the well-known propagator for µ = 2 reads where K α ≡ K 2 α . In result (13), H 1,0 1,1 [·] stands for a Fox H-function [87]. For 0 < α < 1, the above propagator displays a non-differentiable peak (cusp) at the origin. The solution (13) is equivalent to the series representation Employing standard theorems for the Fox functions [87] one can show that for |y| ≫ √ K α t α the following asymptotic stretched Gaussian behaviour emerges [23]: From equation (13) (or from symmetry arguments) it is immediately clear that odd moments vanish. In turn, the behaviour of the even moments is strongly influenced by the stretched Gaussian behaviour (15). One finds [23] Higher order integer moments can be expressed in terms of the variance y 2 as follows, For α = 1 (normal diffusion) one recovers the moments characterising the typical Gaussian propagator, In particular, equation (17) can be used to calculate the kurtosis. This quantity is a measure of the "tailedness" of a given probability distribution, defined as Recall that for a normal distribution one has β 2 = 3 in one dimension. In contrast, a fat-tailed distribution exhibits a large skewness or kurtosis, relative to that of a normal distribution. Distributions with β 2 > 3 are called leptokurtic (as opposed to distributions with β 2 < 3, which are termed platykurtic). For the particular case µ = 2 described by result (16) one finds Thus, β 2 decreases monotonically from β 2 = 6 for α = 0 to β 2 = 3 for α = 1 as α increases. In other words, the tails become less fat with increasing α, and in this sense the distribution becomes more Gaussian-like. However, for any value α < 1, the distribution remains leptokurtic with β 2 > 3.

Growing domain
Next we compare the behaviour of positive integer moments up to the kurtosis of the particle distribution for a symmetric CTRW with their counterparts in a uniformly growing domain.
In the growing domain, the coordinate y of a physical point (hereafter also termed "physical coordinate" or "Eulerian coordinate") is no longer stationary, since each physical point is advected by the growing domain. Thus, the distance between a physical point at y and the origin 0 changes in time as the domain expands. In the following, we will also assume that a random walker "sticks" to the physical medium while it does not jump ‡-consequently the walker is also advected by the medium as it expands. It is convenient to describe the time evolution of y in terms of its initial position x ≡ y 0 , hereafter termed "Lagrangian coordinate". One has where a(t) > 0 is the so-called scale factor withȧ > 0 for a growing domain. Note, however, that our formalism also accounts for the case of a shrinking domainȧ < 0. Analogously to section 2.1 we consider the case of a separable random walk with a (time independent) jump length PDF λ(y) ∼ σ −µ |y| −1−µ for |y| ≫ σ and a waiting time PDF ϕ(t) ∼ (t/τ ) −1−α for t ≫ τ . A description of the random walk in terms of the Lagrangian coordinate is especially well-suited, since the kinetics can then be mapped onto the original domain, at the expense of having to deal with a time-dependent jump length PDF. Indeed, from probability conservation, the jump length PDF on the fixed domain (probability per unit length to take a jump of length x) then reads λ(x|t) = a(t)λ(y = a(t)x).
Assuming that the Fourier transform λ(k) ≡ F [λ(y)] takes the asymptotic form (5) for small wave numbers k, one finds that the Fourier transform λ(k x |t) ≡ F [λ(x|t)] with respect to the Lagrangian coordinate x is given by (22) ‡ For instance, a tracer bead intermittently stuck in an expanding hydrogel.
Using both equations (22) and the asymptotic form (4) of the Laplace-transformed waiting time PDF, a formalism similar to the one introduced in section 2.1 leads to the corresponding description in terms of an FDE for the evolution of the particle's PDF W 0 (x, t) on the original domain (for details we refer to [59]), Note that the main difference with respect to the case of a static domain given by result (8) is that the anomalous diffusion coefficient is now multiplied by the time dependent prefactor a −µ (t). In other words, one has an effective, time dependent anomalous In the case of a growing (shrinking) domain, one can therefore interpret that the diffusive steps measured in Lagrangian coordinates become shorter (larger). In general, this time-dependent diffusion coefficient complicates the solution of equation (23) considerably. When α = 1, an analytic solution can be obtained for any value 0 < µ ≤ 2, which includes the parameter range 0 < µ < 2 describing Lévy flights (see section 4). In contrast, a solution in closed form does not appear to exist when α = 1. However, for µ = 2 a careful analysis of the moments of the distribution suffices to unveil a drastic change in the behaviour of the solution with respect to the case of a static domain. For this particular case, general expressions for the Lagrangian moments x n 0 are available, whence expressions for the Eulerian moments y n 0 immediately follow via the relation y n 0 = a n (t) x n 0 .
2.3. Behaviour of the moments for 0 < α ≤ 1 and µ = 2 Our starting point is equation (23), which in the present case becomes with K α ≡ K 2 α . Moments of different order can be obtained by multiplication with the spatial variable x raised to the corresponding power and by subsequent integration.

Variance
For a generic scale factor a(t), one has the formula [59] In the case of a power-law expansion a(t) = (1 + t/t 0 ) γ (with γ ≥ 0, whereby γ = 0 corresponds to the case of a non-growing domain), one has where 2F1 (·) denotes the regularised hypergeometric function. From relation (15.3.7) on page 559 of reference [88] one directly finds when z → ∞. Similarly, using relation (15.3.13) on page 560 of the same reference one obtains Equation (27) in combination with equations (28) allows one to conclude that the Lagrangian variance displays three different asymptotic long-time regimes, depending on the specific values of α and γ, namely [59] for γ < α/2, for γ = α/2, and for γ > α/2. In this latter case, x 2 0 tends to a constant value as t → ∞, implying that the Lagrangian propagator "freezes" as a result of the fast domain growth.
Correspondingly, in physical space one gets for γ < α/2, for γ = α/2, and In the first case γ < α/2, the domain growth is slow enough to ensure that the long-time behaviour of the variance is the same as in case of a static domain, except for the fact that one has a modified effective diffusion coefficient K α α/(α − 2γ). In contrast, for γ > α/2, the particle drift associated with the deterministic domain growth (the so-called "Hubble drift" in the language of cosmology) becomes so large that the intrinsic diffusive motion of the particle only represents a negligible perturbation, and thus y 2 0 ∝ t 2γ at sufficiently long times. Finally, in the marginal case γ = α/2 the asymptotic variance displays the same qualitative time dependence as for a non-growing domain, but a logarithmic correction appears.
Since a sufficiently fast power-law growth γ > α/2 implies that the long-time behaviour is essentially dominated by the Hubble drift, this will obviously also hold true for faster types of domain growth. An interesting example is the case of an exponential growth a(t) = exp(Ht) with H > 0. This yields implying y 2 (t) 0 ∝ exp(2Ht), as expected.

Fourth-order moment
In contrast to the case of a static domain, on a growing domain the fourth-order moment is related to the variance in a more intricate fashion [59], Using formula (26) and performing the corresponding fractional derivative, one obtains For the specific case of a power-law growth, one has No elementary analytic expression for the above integral appears to exist. However, it can be evaluated numerically. Conversely, equations (28a) and (28b) allow one to infer the long-time behaviour of the fourth-order moment. As in the case of the second-order moment, three different regimes can be distinguished. For a slow expansion γ < α/2, one has In the marginal case γ = α/2, we find Finally, when γ > α/2, since x 4 0 ∝ t/t 0 0 dzz α−2γ−1 , in the long-time limit x 4 0 tends to a constant value (to be determined numerically). As was the case for the variance, the Lagrangian fourth-order moment tends to a constant value. In fact, it can be proven that whenever a "freeze-out" of the variance takes place, it propagates to all even-order moments (to this end, one can e.g. take equations (62) and (63) in reference [59] as a starting point). Of course, this can only mean that the Lagrangian propagator tends to a stationary profile as t → ∞.
For the sake of completeness, we also give here the result for an exponential growth a(t) = e Ht with H > 0. One finds where 0 F 1 (·) stands for the confluent hypergeometric function. The integral on the right hand side remains well-defined in the limit t → ∞ and the asymptotic value of the fourth-order moment is

Kurtosis
Let us briefly discuss the behaviour of the kurtosis (19) on the basis of the results for the fourth-order moment. To start with, note that this quantity refers to the form of the propagator, and is therefore scale invariant as long as the domain growth is uniform (the only case we consider throughout the present work). Thus, i.e., one may indistinctively use Eulerian or Lagrangian coordinates for the computation of the kurtosis. For a symmetric walk, this gives The main novelty introduced by the domain growth is a time-dependent kurtosis in the subdiffusive case 0 < α < 1 (in the normal diffusive case α = 1, the kurtosis remains a stationary quantity). This implies that the propagator for the present case cannot be obtained by a simple rescaling of the propagator referring to a static domain. More precisely, the difference β static 2 − β 2 grows in time and attains a maximum value β static 2 − β ∞ 2 as t → ∞. For the special case of the power-law growth studied previously, the behaviour of the hypergeometric functions in the case γ ≤ α/2 yields (cf. equations (29a) and (35a)) In particular, this implies β ∞ 2 > 3 for γ < α/2 and β ∞ 2 = 3 for γ = α/2. Starting from the case γ = 0 of a non-growing domain and increasing γ, as one approaches γ = α/2 from below the final value of kurtosis decreases monotonically, and the form of the final propagator increasingly resembles a Gaussian. In fact, when γ = α/2 the stationary kurtosis takes the value 3, i.e., that of a Gaussian PDF. In this sense, even though the solution of equation (25) remains non-Gaussian in this case, one may still speak of Gaussian-like behaviour.
The physical origin of the time decrease of the kurtosis observed in the case α < 1 with 0 < γ ≤ α/2 (and also for γ > α/2, see below) is intriguing; for α = 1, equation (25) describes scaled Brownian motion, and the propagator remains strictly Gaussian at all times. In contrast, for α < 1 the propagator is fat-tailed at short times; however, as time goes by, the effect of the Hubble drift on these fat tails of the distribution appears to become stronger than in the central part. This could explain that for α < 1 the kurtosis takes values which are increasingly close to the Gaussian value β 2 = 3 in the course of time. We have no physical explanation for the Gaussian-like behaviour observed when γ = α/2, other than the fact that this particular value separates the diffusion-dominated regime from the regime dominated by the domain growth [59].
In figure 1 we show simulations results for a power-law expansion a(t) = (1 + t/t 0 ) γ with t 0 = 10 3 and different values of γ ≤ α/2. The curves for β 2 (t) can be seen to approach the value of β ∞ 2 given by result (40). In particular, the decrease of β ∞ 2 with increasing γ predicted by equation (40) is confirmed. The theoretical results compare very favourably with Monte-Carlo simulations at times sufficiently long to reach the diffusive regime.
In figure 2 we illustrate how the form of the propagator changes as its kurtosis evolves in time. We show the Lagrangian propagator for a power-law expansion a(t) = (1 + t/t 0 ) γ with t 0 = 10 3 and γ = α/2 = 1/4 at two different times. In panel (a), we display the propagator at a comparatively short time (t = 2 14 ). We have been able to solve the FDE computationally up to this time (details of the algorithm can be found in reference [89]). The propagator is still quite pointy at x = 0 and thus not too different from the shape on a static domain (γ = 0). In more quantitative terms, at t = 2 14 (25) obtained via a fractional finite-difference method [89] with spatial discretisation ∆x = 0.1 and time discretisation ∆t = 0.1 (solid line). The dashed line represents the exact solution for the static case at t = 2 14 . In panel (b) we show the propagator at t = 2 34 obtained from simulations (empty circles), and we compare it with a Gaussian whose variance is taken to be that of the real distribution.  Figure 3. Kurtosis for a subdiffusive random walk with α = 1/2 and K α = 1/2 on fast growing domains (γ > α/2 for power-law expansions and H > 0 for exponential expansions). The two upper curves correspond to an exponential scale factor a(t) = exp(Ht), whose logarithmic derivative is, respectively, H = 10 −5 and 10 −6 . The two bottom curves correspond to a power-law scale factor a(t) = (1 + t/t 0 ) γ with t 0 = 10 4 and γ = 3/2 and 3/4, from top to bottom at t = 10 8 . Symbols represent simulation results from 10 6 realisations. Numerical solutions are depicted by solid lines. The value of β 2 (t) was computed by numerical integration of the Lagrangian fourthorder moment and subsequent division by the analytical expression of the squared Lagrangian second-order moment. The dashed line depicts the value of β static 2 ≃ 4.71 (cf. equation (20)).
In contrast to the case γ ≤ α/2 for a sufficiently fast domain growth (γ > α/2) the final value of the kurtosis β ∞ 2 increases as γ grows. As shown in reference [59], in this γ-regime the Lagrangian propagator eventually freezes and its final form is always leptokurtic (β ∞ 2 > 3). This also implies the long-time freeze-out of all the associated moments, which tend to finite values as t → ∞. In particular, this holds for the secondand the fourth-order moments, from which the kurtosis is computed. As explained in reference [59], the physical reason for the observed freeze-out is the irrelevance of the diffusive spreading with respect to the Hubble drift at sufficiently long times. As a result of this, after a characteristic time t char > t 0 (see figure 3), changes in the form of the Lagrangian propagator and in the associated kurtosis become negligibly small, and the monotonic time decrease of the kurtosis saturates at a value β ∞ 2 ∈ (3, β static 2 ). Of course, t char will depend on both t 0 and γ. For larger values of γ, one expects a decrease of t char , since the Hubble drift becomes dominant with respect to the diffusive spreading at earlier times, and hence the saturation in Lagrangian space becomes faster. This entails a stronger memory of the tailedness displayed by the early-time propagator, and therefore a larger γ leads to a larger β ∞ 2 . Note that this is just the opposite of what happens when 0 < γ ≤ α/2. Our findings for the case γ > α/2 are fully confirmed by the results shown in figure 3, which displays a comparison between theory and simulations in the long-time regime.
Summarising, Gaussian-like behaviour is only observed when γ = α/2. For any other value the kurtosis of the final distribution is always > 3, i.e., the final PDF is leptokurtic and a strong signature of the early-time PDF persists for arbitrarily long times. Note, however, that the asymptotic value β ∞ 2 does not depend on the characteristic time t 0 , which only has an influence on the transient behaviour.
Of course, in the exponential case a(t) = e Ht with H > 0, a freeze-out of the Lagrangian propagator also takes place. It is, however, striking that β ∞ 2 = 3 × 2 1−α , regardless of the value of H. For any α < 1 one has β ∞ 2 > β Gaussian 2 ≡ 3, but β ∞ 2 < β static 2 ≡ 3Γ(α)Γ(1 + α)/Γ(2α)-see equation (20). Results for two different values of H > 0 are displayed in figure 3. The kurtosis β 2 (t) in this case is seen to approach β ∞ 2 at a time of the order of 1/H. As already anticipated β ∞ 2 does not depend on H. Of interest is also the behaviour in the case of exponential contraction, i.e., a(t) = e Ht with H < 0. In this case it is easier to work in physical coordinates. In the limit t → ∞ one has the asymptotic behaviour whence For further details of the calculations of the moments, including the use of Tauberian theorems, we refer to reference [59]. As expected, when α = 1 the kurtosis is equal to three, whereas for α < 1 the kurtosis always grows in time (in this case, it is proportional to t 1−α ).
Other types of contraction can also be considered. The case of power-law contraction, i.e., a(t) = (1 + t/t 0 ) γ with γ < 0, is also covered by our formalism. The kurtosis β 2 (t) increases in time, until a limiting value β ∞ 2 is attained. This limiting value is still given by equation (40), which also holds for γ < 0. For a given α < 1, the value of β ∞ 2 grows with increasing |γ|. We close this subsection by noting that the observed time increase of β 2 for shrinking domains is likely to be related to the inversion of the direction of the Hubble drift with respect to the case of a growing domain (in uniformly shrinking domains, the Hubble drift tends to bring any two physical points closer to each other, whereas it tends to separate them in a uniformly growing domain). We actually conjecture that beyond the two cases with exponential and power-law scale factor studied here, the kurtosis of the propagator generated by subdiffusive walks with α < 1 displays a time decrease (increase) on any uniformly growing (shrinking) infinite domain.

Lévy flights
As already anticipated, in the case of Lévy flights (α = 1 and µ < 2) the propagator W 0 (x, t) associated with equation (23) can be explicitly obtained. One has where is a symmetric Lévy-stable density with exponent µ and scale factor σ L . Note that for µ = 2, this Lévy density becomes a Gaussian PDF with standard deviation σ = √ 2σ L . By definition, the second-order moment of a Lévy flight is infinite. However, one can define a typical width as w µ (t) ≡ C µ σ L (t), where C µ is a constant chosen in such a way that holds, where P is a predetermined probability. The typical width w µ will be seen to play an important role when addressing the kinetics of mixing of two initially localised Lévy pulses evolving on a growing domain (cf. section 4).

Fractional diffusion equation in physical coordinates
It is possible to obtain the PDF W * 0 (y, t) to find a particle inside the interval [y, y + dy] at time t by noting that it is related to the PDF W (x, t) for finding the particle in the corresponding interval [x, x + dx], where x = y/a(t). One has [59] W * 0 (y, t) = W 0 (x = y/a, t) a(t) .
This implies the three relations and Inserting relations (47) to (49) into equation (23) one eventually finds the sought equation

t)x, t)] x→y/a(t) .(50)
In the case µ = 2 it is possible to directly obtain the behaviour of the physical moments y n (t) ≡ ∞ −∞ y n W * 0 (y, t)dy by multiplication of equation (50) with y n and subsequent integration over y. This yields a hierarchy of differential equations for the physical moments.

CTRW in a velocity field
The diffusion of a particle in a uniform velocity field can be regarded as the motion of a random walker dragged by a fluid flowing with velocity v with respect the laboratory reference frame S L (for simplicity, the velocity v will hereafter be assumed to be stationary unless otherwise specified). An example would be a frog performing random jumps with statistically distributed waiting times on a wooden log that is longitudinally floating downstream on a river. On a more microscopic scale one could imagine a tracer particle subdiffusing in a hydrogel that itself is slowly streamed in a fluidic device.
Let us now introduce a second reference frame S 0 in which the deterministic contribution of the velocity field is subtracted from the overall particle motion, i.e., a frame which follows the fluid that drags the particle along. Clearly, S L moves with velocity − v with respect to S 0 . Let us respectively denote by W 0 ( y, t) and W ( y, t) the walker's PDF in S 0 and S L . On a static domain the relation between both PDFs will be given by the Galilean transformation W ( y, t) = W 0 ( y − vt, t). In particular one has W (y, t) = W 0 (y − vt, t) in one dimension. In Fourier-Laplace space equation (51) becomes In the case of a growing domain the particle is not only advected by the velocity field but also experiences an additional drift as it is dragged by the physical medium (i.e., the aforementioned Hubble drift). Recalling the previous example of a hopping frog on a log floating downstream, one might wonder what the equation of motion of the frog would be if the log were replaced by a linear rubber strip that expands uniformly with a certain scale factor. Before providing the answer to this question, we will first address the simpler case of a static domain, as done in section 2.

Static domain
From equations (52) and (6) one finds or, equivalently, whence the Galilei-invariant (GI) advection-diffusion equation follows. The operator (∂ t + v ∂ y ) 1−α is the fractional material derivative introduced by Sokolov and Metzler [90]. It is defined by From this expression and from relation (9) one can see that the fractional material derivative reduces to the GL derivative if v = 0. As mentioned in the Introduction, the GI fractional advection-diffusion equation (55) for µ = 2 was recently obtained by Cairoli et al. [85]. While the standard material derivative, corresponding to the limit α = 1 reflects the GI of a standard physical system, the power (1−α) reflects the spatiotemporal coupling in a waiting time-random walk scenario with a constant relocation speed v.

The fractional derivative in direct (position-time) space is [91]
∂ ∂t + v ∂ ∂y Here, it is assumed that W (k, t) satisfies the condition (11) of good behaviour at t = 0. § The fractional material derivative can be rewritten in terms of the RL integral in the convenient form (57) is slightly different from that in [91], since in that reference the Fourier transform is defined as in equation (3) but with the replacement k → −k).

In this way the GI advection-diffusion equation (55) becomes
Note that in the above expression the time fractional derivative is applied to W 0 , which corresponds to the particle distribution in the reference frame S 0 . This frame moves by a distance vt during the time t. As a result of this the evaluation point for the derivative is shifted by the same quantity. Finally, it should be noted that the propagator (53) and the corresponding FDAE (59) differ from those considered in references [92] and [93]. In particular, for µ = 2, the propagator studied in these works can reach non-physical negative values with moments that are only correct up to order two.

Growing domain
We continue to assume that the particle is subject to the influence of a constant velocity field (the velocity measured with respect to the laboratory frame S L is v). As in the case of a symmetric walk it is convenient to work in Lagrangian coordinates. If the domain expands uniformly with the scale factor a(t), the Lagrangian distance travelled by the reference frame S 0 (the frame moving with velocity v with respect to S L ) during the time t is Here, the transformed time T (t) replaces the physical time t as a result of the difference between the Lagrangian distance Λ(t) travelled by S 0 and the distance vt that it would cover on a static domain. For instance, an exponential growth a(t) = exp(Ht) with H > 0 yields, whereas a power-law growth a(t) = (1 + t/t 0 ) γ gives and Thus, for a sufficiently fast growth (exponential or power-law with γ > 1) one has an asymptotic finite value T ∞ ≡ T (t → ∞), and correspondingly the asymptotic Lagrangian distance Λ ∞ ≡ Λ(t → ∞) is also finite.
Notice that in cosmology, Λ, T and ν are known as "comoving distance", "conformal time", and "comoving velocity", respectively [94]. The PDF W (x, t) in the laboratory frame S L is just the PDF W 0 (z, t) in the comoving reference frame S 0 , but with the shifted position z = x − Λ(t). In other words, In figure 4 the above relation is illustrated for a stretching domain with a(t) = (1 + t/10 3 ) 1/8 by comparing simulations results for W 0 and W . The Lagrangian propagators are also compared with the corresponding ones for the case of a static domain. Note that the width of the propagator is decreased with respect to the static case, since the jump length measured in Lagrangian coordinates is divided by a(t). Now, as we already know from section 2, W 0 (x, t) satisfies relation (23). Since Λ(t) may be a complicated function of t, in general no relationship between the Fourier-Laplace transforms of W (x, t) and W 0 (x, t) similar to equation (52) can be obtained from result (63). Therefore, we can no longer use the straightforward procedure of section 2 to derive the FDAE that we seek. For this reason, we will follow a different path involving equations (23), (59), and (63).
In view of equations (59) and (23) a reasonable guess for the FDAE is since Λ(t) and vt in equations (64) and (59) are, respectively, the displacement in the laboratory frame S L of the comoving frame S 0 during the time t. We confirm that equation (64) is indeed the correct FDAE by showing that W (x, t), as given by equation (63), satisfies relation (64). First, let us evaluate the left-hand side of equation (64), . (65) Next, we evaluate the right-hand side of equation (64), where, in the last step, equation (23) was taken into account. Comparing equation (65) with (66) we conclude that W (x, t) = W 0 (x − Λ(t), t) indeed satisfies relation (64). Defining the FDAE (59) for a uniformly growing domain can be rewritten in a way similar to equation (55), Equation (68), or equivalently equation (64) (see just below), is one of the main results of this paper. Comparing equations (64) and (68) one can see that both are equivalent if In order to prove this we first note that equation (67) can be rewritten as (see equation (12)) Conversely, for any differentiable function G(x, t), one has Taking G(x, t) = 0 D −α t W (x+Λ(t), t) and assuming that (d/dt) 0 D −α (see equation (10)), one finds that the right hand side of equations (69) and (70) coincide, implying that relation (69) indeed holds. Even though the velocity v has been assumed to be constant, it is worth noting that equations (64) and (68) still hold when the reference frame S 0 moves with a timedependent velocity with respect to S L . In this case, one simply replaces ν = v/a(t) with ν = v(t)/a(t).
We close this section by deriving the FDAE in terms of the physical coordinate y. The derivation proceeds along the same lines as in the field-free case (cf. section 2). The relations (46) to (49) for this latter case are completely analogous to the equations relating W * (y, t) and W (x, t) in the presence of the velocity field. Indeed, one has and, correspondingly, as well as and Inserting equations (72) to (73c) into (64) and taking into account equation (63) one eventually obtains the result where Let us once more recall that equation (74) holds for well-behaved functions W * (y, t).
Strictly speaking, the RL-fractional derivative therein must be replaced with the GLfractional derivative 0 D 1−α t , and consequently, 0 P 1−α t must also be replaced with the corresponding operator, defined via the equation

Moments of the propagator
From the exact relationship (63) one can easily find the moments x n of W in terms of the moments x n 0 of W 0 . One has Restricting ourselves to the case µ = 2 and 0 < α ≤ 1 (subdiffusive case), the secondand fourth-order moments of the PDF are respectively given by equations (26) and (33).
By symmetry, odd moments vanish trivially, x 2n+1 0 ≡ 0. Taking all this into account in equation (77) we find explicit expressions for the first four Lagrangian moments: Focusing on the first two equations, equation (78a) tells us that the first moment x is simply given by the deterministic shift, i.e., the Lagrangian distance Λ(t) travelled by S 0 after a time t. Conversely, equation (78b) implies that x 2 0 − x 2 0 equals x 2 − x 2 , i.e., the variance in S 0 is the same as in S L . Thus, regardless of the value of α the dispersion properties of the random walk are not affected by the velocity field. In contrast, in the case of a constant external force [61] this only holds when the particles are Brownian (α = 1). As soon as 0 < α < 1, the dispersion properties are altered both on a static domain [23] and on a growing domain [61].
It is instructive to check that the expressions for the moments obtained above can be directly recovered from equation (68). In order to deduce equations (78a) to (78c) from equation (68) we first multiply this latter equation with x n and subsequently integrate over space. We are then left with the hierarchy where Performing the change of variable z = x−Λ(t)+Λ(τ ) in this integral, using the binomial expansion for z n , and integrating by parts, one finally obtains where For n = 1 equation (81) becomes d x /dt = ν, and equation (78a) follows immediately. For n = 2 relation (81) becomes 2,2 .
From this equation and from result (78a) one immediately finds equation (78b) by taking into account that dJ We now proceed to derive equation (78c) from result (81). For n = 3, relation (81) can be written as after a straightforward calculation making use of equation (78a). Inserting equation (78b) into (84) one obtains Equation (78c) then follows immediately from (85). The validity of result (78d) for the fourth-order moment can also be corroborated in a similar way. The calculation is straightforward but not shown here explicitly.

Mixing of diffusive pulses
We now proceed to study the influence of a velocity field on the mixing properties of two diffusive pulses on a uniformly growing one-dimensional domain. We consider both cases of normal and anomalous diffusion. For convenience, our analysis will be carried out in Lagrangian coordinates.
As mentioned in the Introduction, mixing properties are crucial to understand encounter-controlled reactions involving pairwise interactions, as originally modelled by Smoluchowski [95]. In biological cells, for instance, monomers of regulatory proteins meet diffusively and form dimers [96], and the dimer then diffuses to its designated binding site on the genome or a DNA plasmid [97,98]. ¶ Analogous processes need to run off in vesicles designed as artificial cells [99]. In the case of particles advected by the growing domain, for instance, by the expanding cytoskeleton in a living biological cell or compartments in growing vesicles [99], it is intuitively clear that the associated Hubble drift will reduce diffusional mixing and thereby decrease the reaction rate [64][65][66]. For sufficiently fast domain growth the lack of mixing stems from a premature freezing of the pulse propagators [63][64][65][66], since the Lagrangian step lengths become increasingly short in the course of time (for a contracting domain, the Lagrangian steps become larger and one has the opposite effect).
In our problem, the random motion of each pulse is not only subject to a Hubble drift, but also to an additional drift arising from the velocity field. In order to formulate the problem in the most general form, we will assume that each diffusing particle "feels" a different velocity field. In other words, the force underlying this velocity field can be thought of as being able to discriminate each particle by a distinctive property. For instance, if the physical origin of the force is an electric field acting on charged, overdamped diffusive particles, then this property will be the sign and absolute value of their electric charge. In the particular case where both particles have the same charge, they will experience the same biasing force, and the mixing problem will be equivalent to its zero-field counterpart, save for a coordinate shift. In the biological context, we could be thinking of growing neuron cells, in which messenger RNA molecules are shuttled along by molecular motors [83]. Depending on the orientation of the molecular track these motors are walking, their direction may be towards either extremity of the pseudoone-dimensional cell.
More specifically, consider two walkers labelled with indices 1 and 2. Let x (1) 0 and x (2) 0 denote their initial positions. Without loss of generality we assume that x (1) 0 < x (2) 0 . Owing to the external force acting on each walker, the maxima are shifted with respect to their initial positions as the pulses associated with each walker widen. Correspondingly, one has x (1) The velocities v (1) and v (2) will hereafter be considered to be constant for the sake of simplicity.
The two-pulse PDF can be written as a normalised linear combination of the propagators for the respective initial conditions, i.e., Here, one encounters the difficulty that if at least one of the particles is subdiffusive, its propagator is not known, and hence its contribution to the joint PDF is also unknown. Nevertheless, it is possible to carry out a semiquantitative study of particle mixing on the basis of the second-order moments. To this end, let us first introduce the Lagrangian half-width of a zero-field single-particle propagator as Let us further define two characteristic points x (1) M + w (1) (t) and x (2) (2) (t) where, following the notation of equation (87), w (1,2) (t) > 0 denotes the Lagrangian half-width of the symmetric propagator W (1,2) 0 (x, t). Mixing after a time t will be considered to be weak (in a statistical sense) if x In contrast, when we will simply speak about "mixing". Similarly, if "<" can be replaced with "≪" in equation (89) we will speak about "strong mixing". The last two terms on the righthand side are positive and independent of the respective velocities. Therefore, for a given set of diffusive properties, the mixing behaviour depends on the first term, which is influenced by the relative velocity v (1) − v (2) and by the scale factor a(t) entering the definition of T (t). Clearly, a positive (negative) relative velocity favours (hinders) mixing, as is the case in the particular case of a static domain (this latter case is recovered by setting T (t) ≡ t).
Note, however, that the domain growth introduces a key difference in the behaviour. On a static domain the strong mixing condition is always fulfilled provided that one waits long enough, irrespective of whether the diffusive pulses are normal or anomalous (mixing becomes maximal when both pulse peaks overlap, i.e., when the distance d M ≡ x (2) M − x (1) M between the two peaks vanishes). However, for a given initial pulse separation and relative velocity, it is intuitively clear that a sufficiently fast domain growth may freeze both pulses before their mixing becomes significant. More precisely, for v (1) > 0 and v (2) < 0, one has weak mixing at arbitrarily long times if where T ∞ ≡ lim t→∞ T (t) < ∞ and w (1,2) ∞ = lim t→∞ w (1,2) (t) < ∞ are the long-time asymptotic values of the half-widths. Thus, a sufficiently fast domain growth favours the localisation of the Lagrangian propagators about the respective initial conditions, thereby preventing that the memory of the latter is eventually lost by mixing.
In terms of physical coordinates, the situation described by equation (90) corresponds to the case when the Hubble drift is so strong that the pulses separate from each other at a rate much faster than the typical growth rate of their half-widths. Therefore, the overlap of both pulses remains negligibly small at all times.
In what follows, we focus on the case of a set of two particles with identical diffusive properties (w ≡ w (1) = w (2) ) but subject to opposite drift velocities v ≡ v (1) = −v (2) . Without loss of generality the midpoint between the two pulse peaks will be chosen as the origin, i.e., x 0 ≡ x
In order to study the mixing kinetics the parameter P in equation (45) should be chosen large enough to ensure significant mixing as soon as both tails overlap. In general, for a given value of P the corresponding value of C µ (as well as the associated characteristic width w µ (t) = C µ σ L (t)) must be computed numerically. However, for the Gaussian case µ = 2 one can obtain an analytic expression, namely, C 2 = 2erf −1 (P ). This result is in agreement with what we had already anticipated for the Brownian case since, for a Gaussian distribution, the probability that a particle is found within an interval of half-width w(t) = 2σ(t) = 2 √ 2σ L (t) is P = 0.9545, i.e., precisely the value which follows from the relation erf −1 (P ) = √ 2. In the case of two Gaussian pulses with identical diffusive properties and opposite drift velocities this last result implies that when the weak mixing condition (88) holds the overlap of both tails as given by remains below 5% at all times. For a power-law scale factor a(t) ∼ t γ the different possible subcases are given in Table 1. From this table one concludes that for v = 0 and x 0 ≪ w ∞ strong mixing const. const. Table 1. Asymptotic long-time behaviour of T (t) and w(t) for a power-law scale factor with characteristic exponent γ. The results for T (t) and w(t) stem, respectively from the long-time behaviours of equations (62) and (87).
does not occur when γ > α/2. However, a nonzero value of v may completely change this scenario and bring about strong mixing for sufficiently long times. For v > 0, when α/2 < γ ≤ 1, one has T ∞ = ∞, implying that the crossing of the two maxima will eventually occur with certainty. In contrast, when the domain growth is fast enough to ensure that T ∞ < ∞ the propagator evolves towards a steady state. In this latter case strong mixing will never take place if it has not already occurred by the characteristic time t C at which W (x, t C ) can be considered to be practically indistinguishable from W (x, ∞). According to equation (89) strong mixing may be observed at a finite time t if the initial separation distance is small enough-more precisely, if x 0 ≪ vT (t) + w(t) holds. Figure 5 displays the temporal evolution of two subdiffusive pulses (α = 1/2) on a growing domain with scale factor a(t) = (1 + t/10 3 ) 3/4 . Since γ > α/2 we can see that the probability of overlap grows in time. However, the respective pulse widths remain almost stationary throughout the time window spanned by the represented set of plots (the length of this window is ∼ 10t 0 ). Thus, we conclude that in the present case mixing is driven by the opposed velocity fields rather than by the spreading of the pulses. The maximum overlap probability is attained after a time t M , corresponding to the crossing of both pulses, i.e., to a vanishing peak separation d M (t M ) = 0. This characteristic time can be obtained from the solution of the implicit equation x 0 = vT (t M ). Only at this specific time t M does the PDF correspond to the zero field solution for a single particle, For a fixed time t = 10 4 figure 6 displays a series of snapshots of the Lagrangian propagator corresponding to different values of v. For the chosen value of the domain growth exponent (γ = 1/2) and of the anomalous diffusion exponent α the inequality 1/4 = α/2 < γ < 1 holds, and so the mixing is once again driven by the relative velocity term. As in the static case, increasing v results in enhanced mixing. For the chosen parameter set the degree of mixing at a given time is maximised by the velocity value v = 0.162. For this precise value, d M (t = 10 4 ) = 0, i.e., one has a single-hump propagator.
Finally in figure 7 we study a case where the domain growth is so fast (γ = 2) that mixing is absent at all times. As one can see, the pulse widths remain practically the same at all times, and the deterministic displacement of the pulse peaks is greatly slowed down by the domain growth, until both pulses become almost stationary.
The theoretical curves in figures 5-7 (obtained by the numerical integration of the   ) were corroborated by CTRW simulations. As expected, a significant discrepancy occurs at short times, since the number of steps taken by the random walker is not large enough to reach the diffusive limit (see also reference [23]). At longer times the agreement is excellent. γ < min{1/µ, 1} 1/µ = γ < 1 1/µ < γ < 1 γ = 1 const. w µ t 1/µ−γ (ln t) 1/µ const. Table 2. Asymptotic long-time behaviour of T (t) and w µ in the case of a power-law expansion with exponent γ.

Lévy flights
To conclude our study about pulse mixing, let us now focus on the specific case of Lévy flights corresponding to the parameter choice α = 1 and 0 < µ < 2. The solution for a two-pulse initial condition can be easily inferred from the one-particle propagator (43) for the zero-field case via the relation (86). The explicit form of the solution is with [σ L (t)] µ = K µ 1 t 0 a −µ (u)du. Let us once again consider the case of the power-law scale factor a(t) = (1 + t/t 0 ) γ . As it turns out the asymptotic long-time behaviour of the typical width depends on whether 1/µ is larger, equal, or smaller than γ. Table 2 summarises the typical subcases that result from the long-time behaviour of T (t) and w µ (t) for different values of γ and µ.
Thus, for v = 0 the mixing of both pulses can be avoided for arbitrarily long times by choosing γ > 1. However, when v = 0 one must only have γ > 1/µ in order to prevent mixing. The reason is that the respective pulse widths grow as t 1/µ , which is not fast enough to ensure significant diffusive mixing of both pulses when their separation distance increases as t γ (with γ > 1/µ) due to the domain growth.
In contrast, for Lévy flights with µ < 1, the dominant contribution to mixing in the long time limit will stem from Lévy diffusion rather than from the biasing fields. Note that such a regime can never occur when the jump length PDF has a finite variance, i.e., in the normal-diffusive case or in the subdiffusive case.
The most clear-cut situation is found in the range 1 < γ ≤ 1/µ. In this regime, since T ∞ < ∞ the positions of both peaks tend to fixed limiting values. The advection velocity v and the initial pulse separation 2x 0 will determine whether or not such limiting values are attained before the two peaks meet. In either case diffusive spreading remains dominant with respect to the Hubble drift after a sufficiently long time and the pulse widths grow without bound. Consequently, at long times the mixing process proceeds via the widening of the respective pulses. Both pulses eventually merge into a single one, which still continues to widen. This is precisely what is seen in figure  8, which displays the evolution of the theoretical propagator for two Lévy pulses with characteristic exponent µ = 3/4 spreading on a domain whose growth is controlled by a power-law factor with γ = 5/4 and t 0 = 10 3 . The figure also shows simulations results that are in excellent agreement with the theoretical expression (92).
Finally, let us discuss the behaviour when 1 < 1/µ < γ. Here the diffusive spreading remains the main mixing mechanism over a transient period. However, the Hubble drift eventually becomes dominant, and consequently W (x, t) tends to a stationary (frozen) profile in the long-time limit. Depending on the characteristic parameters (γ, t 0 , K µ 1 , x 0 , and v) the distribution W (x, t) will be single-or double-peaked at given time.

Summary and outlook
We investigated several aspects of an anomalous diffusion processes described by a separable CTRW evolving on a uniformly growing one-dimensional domain. First we studied the behaviour of the first moments in the symmetric case, with special focus on the rich phenomenology of the kurtosis. One of the most remarkable features is that this quantity becomes time dependent in the subdiffusive case. As a result of this, an initially non-Gaussian subdiffusive pulse was shown to exhibit a Gaussian-like long-time behaviour for a suitable parameter choice.
We subsequently considered the effect of a velocity field on this CTRW dynamics. We derived the corresponding FDAE (68), which generalises a previous result valid for the special case µ = 2 on a static domain [85]. Indeed, since our bifractional equation holds also for µ < 2, it includes Lévy flights as a particular case. The mathematical form of the FDAE is rather peculiar and not intuitive even in the case of a static domain, given the very straightforward Galilean transformation (51) between the walker's PDF in the lab frame S L and its counterpart in the comoving frame S 0 .
Taking the above results as a starting point we studied the mixing behaviour of a pair of diffusive pulses in a one-dimensional domain whose time evolution is governed by a power-law scale factor. We focused on the case where the pulses are drifting with velocities v and −v as they spread, the spreading of each pulse being characterised by their respective half-widths. In this scenario a sufficiently fast domain growth was found to largely prevent mixing between a pair of normal diffusive walkers or between a pair of subdiffusive walkers. However, a sufficiently large value of v is able to restore mixing. In the superdiffusive case, the behaviour is more complex. For µ ≥ 1 and a sufficiently large initial separation of the pulses mixing is essentially controlled by the relative velocity introduced by the fields, as is the case for normal diffusive or for subdiffusive walkers. In contrast, when µ < 1, diffusional mixing dominates over the deterministic mixing induced by the velocity fields.
The present work may be extended in several directions. One such direction should consider the mixing behaviour of pulses for time dependent velocity fields v = v(t). Another interesting generalisation concerns the case of non-separable CTRWs, e.g., Lévy walks. Finally, one could explore the effect of nonuniform domain growth [62] introducing a spatial dependence in the scale factor. well as the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) for an Alexander von Humboldt Polish Honorary Research Scholarship.