Emergent anomalous transport and non-Gaussianity in a simple mobile–immobile model: the role of advection

We analyse the transport of diffusive particles that switch between mobile and immobile states with finite rates. We focus on the effect of advection on the density functions and mean squared displacements (MSDs). At relevant intermediate time scales we find strong anomalous diffusion with cubic scaling in time of the MSD for high Péclet numbers. The cubic scaling exists for short and long mean residence times in the immobile state τim . For long τim the plateau in the MSD at intermediate times, previously found in the absence of advection, also exists for high Péclet numbers. Initially immobile tracers are subject to the newly observed regime of advection induced subdiffusion for short immobilisations and high Péclet numbers. In the long-time limit the effective advection velocity is reduced compared to advection in the mobile phase. In contrast, the MSD is enhanced by advection. We explore physical mechanisms behind the emerging non-Gaussian density functions and the features of the MSD.


Introduction
One of the simplest equations to describe the transport of tracers in subsurface aquifers (water-bearing layers or permeable rock or sediment) is the advection-diffusion equation [1,2] where G(x, t) denotes the probability density function of a tracer particle, v the constant advection velocity, and D the diffusion constant. The initial condition is G(x, 0) = δ(x), corresponding to a 'point injection' in geoscience. The probability density function encoded by equation (1) with initial condition G(x, 0) is given by the Gaussian with similarity variable x − vt, where the first moment ⟨x(t)⟩ =´∞ −∞ xG(x, t)dx = vt and the mean squared displacement (MSD) ⟨[x(t) − ⟨x(t)⟩] 2 ⟩ = 2Dt are linear at all times. The latter does not depend on v. The alternatively used mobile-immobile model (MIM) is a more elaborate model that takes pores into account, where a tracer can remain immobile for an exponentially random duration [3][4][5]. This linear first order mass transfer is often used to model sorbing solutes, as well [6][7][8]. In an MIM the concentration is split into a mobile density n m (x, t) and an immobile density n im (x, t). In the mobile state tracers are subject to advection and diffusion in the same way as in the advection-diffusion equation (1). MIMs have been used extensively in geophysical systems [3,4,[8][9][10][11][12][13][14][15][16][17][18]. Apart from geophysical applications, we mention two other applications, where motion interrupted by transient immobilisations has been studied. The first application pose biological systems, such as potassium channels in the membrane of living cells [19] and transcription factors [20]. Many systems such as tau proteins [21,22], synaptic vesicles [23], complexes at the endoplasmic reticulum [24], the drug molecule doxorubicin in silica nanoslits [25] and DNA binding proteins [26][27][28][29][30][31][32][33][34][35] may be described in terms of an MIM, in which the residence times in the immobile state are distributed exponentially 5 . The corresponding experiments were conducted in flow cells [26,[32][33][34][35] or occur in live biological cells [27][28][29][30]. In such cellular systems advection may arise due to the action of molecular motors causing streaming in the cytoplasm [36]. Further examples of systems with advection, where in addition immobilisations occur, include DNA molecules in microfluidic setups and bio sensors [37,38]. The second application, featuring transient immobilisations concerns charge carriers in semiconductors, where a recent focus lies on exciton diffusion in layered perovskites and transition metal dichalcogenides [39][40][41][42]. Often, MIMs are not formulated in terms of mean residence times but with a single rate for mass exchange and a solute capacity coefficient that takes different volumes of the mobile and immobile volume into account [3,7,11,12,43]. The moments for mobile, immobile and total density of the MIM have been calculated for various initial conditions while including effects of advection and diffusion [6,7,13]. We here focus on the densities and MSDs at relevant intermediate time scales and unveil interesting new properties in the transport dynamics. In the short time limit initially mobile tracers behave like free Brownian tracers (2) with an MSD ∼ 2Dt. At long times the MSD is linear, with the effective diffusivity D eff containing a term proportional to v 2 , that can yield D eff > D [7,13]. This is in contrast to the solution (2) of the advection-dispersion equation (1), in which the MSD does not depend on v. Below we provide a physical explanation for how this enhanced effective diffusion coefficient is brought about. Another model often used to describe the motion of tracers with immobile periods is the continuous time random walk (CTRW) [39,44,45], for which it was shown that exponential tails emerge in the position density [46] when a drift is present [13,18,[47][48][49][50]. CTRWs have also been analysed for systems with two states, characterised by two waiting time distributions, such that a specific distribution is chosen alternatingly or randomly [51,52]. A similar approach to modelling motion interrupted by transient immobilisations was studied in [53]. In the present work we consider a MIM in the formulation with mean mobile residence time τ m and mean immobile residence time τ im similar to our previous work [54] in the absence of a drift. In [54], the mean immobile residence time τ im and mean mobile residence time τ m are well separated, τ im ≫ τ m , corresponding to the one-dimensional motion of tau proteins without advection. At intermediate times τ m ≪ t ≪ τ im a Laplace distribution of positions const × exp(−const|x|)/2 with fixed variance was shown to emerge whose prefactor depends on the initial condition, and the MSD of initially mobile tracers displayed a plateau at intermediate times.
The main goal of the present work is to analyse how the Laplace distribution and the plateau in the MSD change when advection is added. The transition from the Brownian MSD 2 Dt at short times to 2D eff t implies a crossover regime, in which the MSD grows faster than linear given our finding D eff > D. In fact, we find a sustained superdiffusive regime with a cubic anomalous diffusion exponent in the MSD at relevant intermediate time scales. For low advection velocities we recover the model from [54]. Therefore, to highlight new features, we focus on high Péclet numbers Pe = v 2 τ m /(2D) ≫ 1. An application of the MIM with a high Péclet number may occur for sufficiently long times in subsurface aquifers, in the hyporheic zone or in microfluidic setups. A specific example is the motion of biomolecules in bio sensors, as schematically depicted in figure 1. The biomolecules are inserted into a flow cell and bind to the surface reversibly [38]. The surface is coated with receptors that specifically bind to the molecule of interest. Only the bound molecules can be detected using e.g. surface plasmon resonance [38]. In our model we assume the detector to be completely covered with receptors and consider concentrations well below saturation.
In the next section we introduce the model equations and show two ways to solve the model. The direct way using Fourier-Laplace transform yields the densities in Laplace space and exact expressions for the MSD. Additionally, we show how to solve the advection-MIM using a subordination approach, that produces a physical explanation for the additional term in D eff due to the variance of times the tracers spend in the mobile state. As mentioned in the following, we consider strong advection and obtain asymptotic expressions for the density and the MSD in the presence of a clear time scale separation, i.e. τ m ≪ τ im or τ im ≪ τ m . For clarity we focus on the detailed behaviour of the MSD of the total density, while the results for the immobile and mobile fractions are summarised in the appendix. The dimensionless form of the model depends on the ratio τ im /τ m and the Péclet number only. We use these variables in a phase diagram to analyse the anomalous diffusion in the full parameter space including small Péclet numbers. Specifically, in section 2, we formulate and solve our model. In section 3 we consider tracers that are initially mobile and obtain asymptotic expressions of the density functions and MSDs. Special focus is put on finding the parameter regimes where non-Gaussian displacement distributions and anomalous scaling of the MSD emerge. In section 4 we repeat the same analysis for initially immobile tracers. Finally, we summarise and conclude in section 5. The appendix provides details on the calculations and additional figures of the MSDs.

Formulation of the model
We employ the MIM with mean mobile residence time τ m and mean immobile residence time τ im in an (effectively) one-dimensional setting with position variable x, where n m (x, t) and n im (x, t) denote the line densities of mobile and immobile tracers, respectively. Advection with velocity v and dispersion with diffusion constant D is exclusively affecting the mobile density. A single tracer switches between the mobile state and immobile state following a two state continuous time Markov process, i.e. it follows a Poissonian switching. The realisation (3) of the MIM corresponds to the model used in [54] with the residence time distribution in the mobile state γ(τ ) = exp(−τ /τ im )/τ im . We add a drift in the mobile state here and will show that this has significant consequences. Splitting the total density n tot (x, t) = n m (x, t) + n im (x, t) has the advantage that not always n tot (x, t) is measured in experiments. In a biophysics sensor as sketched in [38], solely the immobile tracers can be measured. In contrast, it is often the mobile density that is measured in geological experiments [12,18]. In both cases it is thus essential to model the two densities separately.
Step or delta injection into the mobile domain is common in geological experiments [10,12,[55][56][57]. Our model (3) is very similar to the MIMs used in geoscience, with the difference that there usually the capacity coefficients and porosity, among others, are used instead of the mean residence times [6,7,13]. We consider the initial condition n m (x, 0) = f 0 m δ(x) and n im (x, 0) = f 0 im δ(x) where δ(x) denotes the Dirac-δ distribution. The factors f 0 m and f 0 im denote the fractions of initially mobile and immobile tracers with f 0 m + f 0 im = 1, which effects the normalisation of the total density, ∞ −∞ n tot (x, t)dx = 1. This corresponds to a single-particle picture with no interactions between the tracers, as is the case in model (3). In this work we consider the cases when all tracers are either initially mobile (f 0 m = 1) or when all tracers are immobile (f 0 im = 1). We briefly discuss an equilibrium fraction of initially mobile tracers in appendix K. Other mixed initial conditions are comparatively uninteresting and do not correspond to the experimental initial conditions we have in mind. The structure of this section is as follows. In section 2.1 solutions are derived using Fourier-Laplace transforms. The solution using a subordination approach is shown in section 2.2. In section 2.3 it is demonstrated that the dimensionless formulation of the model (3) model depends on two parameters only. Finally, we discuss the short-and long-time behaviour in sections 2.4 and 2.5.

Solution in Laplace space
as well as (see equation (A.4)) . The fraction f m (t) =´∞ −∞ n m (x, t)dx of free and the fraction f im (t) =´∞ −∞ n im (x, t)dx of immobile tracers are a function of time and the expressions are given in equation (A.8) in appendix A. From these relations we immediately deduce that the total density is normalised, f m (t) + f im (t) = 1. In the long-time limit, the equilibrium fraction f eq m = τ m /(τ m + τ im ) of all tracers are mobile. We calculate the pth moment (p ∈ N) using where j stands for m, im, or tot [54]. To shorten the notation, we use ⟨x 2 (t)⟩ = ⟨x 2 (t)⟩ tot in the remainder of this work. The lengthy exact expressions for the MSD are given in appendix B. We will study their detailed behaviour below and in sections 3 and 4 for specific initial conditions. The first moment ⟨x(t)⟩ is related to the second moment ⟨x 2 (t)⟩ 0 in the advection-free setting with the second Einstein relation [58] This relation becomes obvious when we look at the moments in appendix B. Since ⟨x 2 (t)⟩ 0 was discussed in [54] in detail, we focus on the MSD here.

Subordination approach
The concept of subordination was originally introduced by Bochner [59] and refers to a process X[τ (t)] with the operational time τ (in many random walk contexts the number of jumps [44]), which has random non-negative increments. For the operational time τ the propagator is known, in our case the Gaussian G(x, τ ) (2). Then, the subordinator relates the stochastic increases of τ to the process (laboratory) time t measured in the real-world observation. In our model the stochasticity comes from the immobilisations, where t increases even while τ is stalling. Assume that we know of a single tracer when it is mobile and when it is immobile. Let i(t) be the index function that is unity if the tracer is mobile at time t and zero otherwise, as shown in figure 2. The index function follows a two-state continuous-time Markov process (telegraph process) with mean residence time τ m in state 1 and mean residence time τ im in state 0. As described in [60][61][62] (compare also [63]), this allows us to write a Langevin equation of the form with the operational time τ , that is a random quantity, and zero-mean white Gaussian noise with correlation ⟨ξ(τ )ξ(τ + ∆τ )⟩ = δ(∆τ ). Figure 2 shows two samples of τ (t) for τ m = 1 and τ im = 2 in the upper panel.
The lower panel shows the corresponding trajectory x(t) for D = 1/2 and v = 2. The propagator with this new variable τ is given by G(x, τ ) (2). In appendix D we show how to obtain the subordinator P(τ, t) and its moments using Laplace transforms. In figure 3 we show the numerical Laplace inversion of P(τ, s) (D.12) for three values of t and the parameters τ m = 1 and τ im = 1/2. In order to verify our approach, we compare the Laplace inversion to the solution P(τ, t) obtained in [23], that is presented in appendix D. We find good agreement between the two approaches, as shown in figure 3. In the long-time limit P(τ, t) converges to a Gaussian with mean µ = tτ m /(τ m + τ im ) and variance σ 2 = 2τ 2 m τ 2 im /(τ m + τ im ) 3 t, as shown by the dashed  line following a Gaussian with that mean and variance (see the moments above). With P(τ, t) we can eliminate τ (t) and find using the well established method of subordination [60][61][62][63][64][65] where G(x, t) denotes the propagator (18) in the mobile phase. Expression (11) can be written in the form (6) of n tot (x, s) by inserting the subordinator P(τ, s) in Laplace space. In appendix D we show how to obtain n m (x, t) and n im (x, t) in a similar way.

Dimensionless form
Our model (3) has the four parameters τ m , τ im , v and D. We now show that in a dimensionless form the model depends only on two free parameters, which significantly reduces the space of parameters we need to consider in order to obtain a full picture of the model. Moreover, this highlights the conceptual simplicity of the model. To this end, we define the new dimensionless variables t ′ = t/τ m and x ′ = x/ √ Dτ m . In these variables the model (3) turns into the set of equations which only depends on the immobilisation ratio τ im /τ m and the Péclet number Pe = v 2 τ m /(2D). The typical length scale for the latter is given by vτ m . The factor 1/2 in the Péclet number is introduced for convenience. To see this, we note that from the Péclet number we obtain the advection time scale τ v = 2D/v 2 , which naturally arises from the solution to the advection-diffusion equation (2) as follows. The typical distance travelled due to advection and dispersion is given by ∆x v = vt and ∆x D = √ 2Dt, respectively. Comparing these distances gives the time scale τ v = 2D/v 2 = τ m /Pe, after which displacements due to advection dominate over diffusive displacements.

Short-time behaviour
In the mobile phase tracers are being propagated with the Gaussian (2). As described above the displacements are diffusion dominated for t ≪ τ v . This means that the tracers follow the same density as in the case without advection, which is described in detail in [54]. To emphasise the effects of advection we therefore choose τ v ≪ τ m , τ im . This defines the short-time limit t ≪ τ v .

Long-time asymptote
The densities of the total, mobile and immobile density are, up to a factor, the same in the long-time limit t ≫ τ m , τ im [5]. In absence of advection, the effective long-time diffusivity is given by [54]. Tracers disperse slower compared to the model of simple diffusion advection without immobilisation (1) with diffusivity D. We now obtain D eff by analysing the asymptotic behaviour t ≫ τ m , τ im of the expressions for the MSD. The exact expressions for the first and second moments are stated in appendix B. These results provide the long-time asymptote of the MSD for all initial conditions, namely, where ∼ denotes asymptotic equivalence with an additional spread ∝ v 2 t compared to the case without advection. Remarkably, the asymptotic dispersion with the new effective diffusion coefficient can hence even be higher than the diffusivity D in the mobile state and overcompensate the slow-down of the spread due to immobile durations. In appendix H we explore the parameter regime that yields D eff > D. As shown in figure H1, the increase of D eff is highest, when τ m and τ im are of the same order. A physical intuition for the additional spread due to advection arises from the special case D = 0, for which x = vτ with the mobile duration τ . Due to the stochastic switching between the mobile and immobile state, an ensemble of tracers will have a distribution of mobile durations τ (the density P(τ, t) in section 2.2), and hence the positions will have a finite spread. This additional spread was derived before [6,13], albeit without a more detailed physical justification. The anomalous transport regime at intermediate times that necessarily has to exist to effect the crossover between the two normal-diffusive regimes, and the important physical consequences will not be explained 6 . We explore the physical mechanism behind the additional dispersion due to advection in the long-time asymptote (14) of the MSD. We start by discretising the advection-diffusion process into discrete steps ∆x i that are normally distributed ∆x i d = N (v∆t, 2D∆t), where each step takes a small duration ∆t ≪ t, τ m , τ im . After n steps the tracer position is given by n i =1 ∆x i = x. In the long-time limit the number of steps n follows a Gaussian n d =N (µ/∆t, σ 2 /(∆t) 2 ), as described in section 2.2 7 . The number of steps n is independent of the steps ∆x i . With the expectation value E and the variance var we obtain the expression [66] var 6 Later we will demonstrate that the transient anomalous diffusion regime can assume substantial time spans with an approximately constant anomalous diffusion exponent. 7 The Gaussian distribution n d = N (µ/∆t, σ 2 /(∆t) 2 ), that only applies in the long-time limit, in principle allows for negative numbers of steps. However, they asymptotically vanish, because the width of the distribution vanishes compared to the mean. We quantify this using the coefficient of variation of step numbers for t ≫ τ m , τ im . Equation (16) demonstrates that in the case of diffusion without advection only the mean mobile duration affects the MSD. In contrast, advection couples to the mean in the first moment and to the variance of the duration that the tracers spend in the mobile state in the MSD (14).
All time scales in the model are finite, and therefore the density converges to the Gaussian (17) in the long-time limit with the effective diffusion coefficient D eff (equation (15)) and the effective advection velocity v eff = vτ m /(τ m + τ im ) (see appendix B.1). In the next two sections we analyse the densities and MSDs for the specific choice of initially mobile and immobile tracers, respectively.

Initially mobile tracers
In this section we assume all tracers to be initially mobile. This situation corresponds to many experimental realisations, when tracers are introduced into the system, such that they have not had a chance to immobilise. We consider the following four aspects. In section 3.1 we assume τ v ≪ τ m ≪ τ im , which corresponds to a high Péclet number and long immobilisations to emphasise the effect of advection compared to [54]. In section 3.2 we consider the density for short immobilisations and high Péclet numbers, τ v ≪ τ im ≪ τ m . The third section 3.3 is concerned with the MSD for which we obtain expressions for superdiffusion at intermediate time scales arising due to advection. In the final subsection appendix I.1 we analyse for which parameters the uncovered superdiffusive regime exists and how it is competing with the plateau regime in the MSD found in [54] for long immobilisations in the advection-free regime.

Densities for long immobilisations
Initially mobile tracers that have not yet immobilised follow a Gaussian distribution corresponding to free Brownian motion without advection for t ≪ τ m . At short times t ≪ τ v advection is negligible in the Gaussian propagator and the density is the same Gaussian with mean vt ≪ √ 2Dt close to zero and variance 2 Dt. At intermediate times τ v ≪ t ≪ τ m advection becomes relevant, and the mobile density takes on the Gaussian form  (5) and (6). We use the parameters τ m = 1, τ im = 100, D = 1/2 and v = 100. These parameters are specifically chosen to be able to resolve the multiple time regimes. In the short to intermediate time regime t ≪ τ m , τ im , where t can be shorter or longer than τ v , the immobile density consists of tracers that immobilised at most once. In equation (F.2) in appendix F we arrive at the asymptotic expression for t ≪ τ m , τ im . The mass corresponding to (19) is´∞ −∞ n im (x, t)dx ∼ t/τ m , corresponding to 10 −3 in the left panel of figure 4. This value is very small and hence the immobile density is not visible in the first panel.
The immobile density in figure 4 appears to be uniform at t = 0.1. Indeed, for τ v ≪ t the short-time density (19) approaches a uniform distribution. Using the properties of the error function we arrive at the asymptotic uniform distribution  (19) for t ≪ τm, where t can be smaller or larger than τ v . In that domain the density consists of the same Gaussian and an additional uniform distribution of immobilised tracers shown by the red area in the second panel. For τm ≪ t ≪ τ im almost all tracers immobilised exactly once and follow a Laplace distribution (23) with different scale parameters for the positive and negative x direction. In the long-time limit t ≫ τ im , after many immobilisations the density follows the Gaussian (17) with mean v eff t and effective diffusion constant D eff . See text for details. Parameters: This shape can indeed be seen in the second panel of figure 4, for which the immobile density remains almost constant for 0 < x < 10. The approximation (19) is shown as the blue dashed line. In the same panel the Gaussian (18) is shown as the dashed red line. The total density is given by The appearance of this regime is new, as compared to the case without advection [54]. A physical picture for the occurrence of the uniform density of immobile tracers is as follows. For τ v ≪ t ≪ τ m advective transport dominates over diffusion. Indeed, the typical distance a tracer moved due to advection is given by the mean position vt, while the typical distance travelled due to diffusion is given by the standard deviation √ 2Dt. The mobile tracers hence move with a narrow Gaussian distribution while a fraction of tracers immobilises with the constant rate 1/τ m , which generates the uniform part of the distribution (see figure 4 for t = 0.1). Now we go to immobilisation dominated intermediate times τ m ≪ t ≪ τ im . This regime is easy to analyse when starting from the density in Laplace space. We have ϕ(s) ∼ 1/τ m for sτ m ≪ 1 ≪ sτ im and find the expression from n tot (x, s) (6), which in time-domain corresponds to for τ m ≪ t ≪ τ im . Expression (23) is a normalised distribution with time-independent parameters, that falls off exponentially in the positive and negative x direction with different coefficients. It is shown in the third panel in figure 4 for t = 10 as the grey dashed line. The density falls of quicker in the direction opposite to the advection velocity, as expected. Noteworthily, the Laplace distribution occurs for all values of the Péclet number, i.e. v can take on any value in the asymptote (23).
For long times t ≫ τ im , we show the Gaussian (17) in figure 4 as a red dashed line for t = 10 4 . The Gaussian (17) contains the effective advection v eff = vτ m /(τ m + τ im ) and the effective diffusivity (15), that we obtained in section 2.5.

Densities for short immobilisations
Now we turn to the case of short mean immobile residence times, τ im ≪ τ m . For short times t ≪ τ v and intermediate times τ v ≪ t ≪ τ im ≪ τ m the mobile density follows the same Gaussian (18) as for long immobilisations. This can be seen in figure 5 in the first panel for t = 0.1. From simulations, we obtain position histograms that we colour-code according to the number of immobilisations. Most tracers have not immobilised at this time, as shown by the dominant black area that follows the Gaussian (18). The tracers that did immobilise follow the same asymptote (19) of n im (x, t) as for long immobilisations, as depicted by the dashed blue line. Beyond τ im immobile tracers become mobile again and contribute to the mobile density. This can be seen in figure 5 in the second panel for t = 1 by the red area in the mobile density denoting tracers with a single immobilisation event. In figure 5, the total density appears to have an exponential tail in the opposite direction of the advection velocity for t = 10 and t = 20 in the third and fourth panel, respectively. We now investigate the origin of this phenomenon. For t ≪ τ m , the fraction t/τ m of mobile tracers is trapped for a short period τ drawn from γ(τ ) = exp(−τ /τ im )/τ im . This is shown as the red area in figure 5 denoting a single immobilisation. Hence, these tracers were mobile for a total period of t − τ . We convolute this with the propagator for advection diffusion. Together with tracers that have not immobilised this gives the total density for τ m ≪ t ≪ τ im , The integral on the right-hand side can be solved and is given bŷ for v 2 τ im ≫ D with the complimentary error function erfc(x) = 1 − erf(x). In appendix G we develop the full expression of the integral (25), that is also valid for v 2 τ im ∼ D, i.e. intermediate Péclet numbers. Approximation (24) is shown in figure 5 as the dashed dark-yellow line and follows the density for x > 70 at t = 10 and matches almost the entire density at t = 20. For x ≪ vt and v 2 τ im ≫ D expression (24) simplifies to This tail approximation is only valid close the origin. For this reason it is not normalised. It is shown in figure 5 as the solid grey line. The tail overlaps with expression (24) shown as the dashed dark-yellow line. The second exponent in expression (26) reveals that the slope of the tail does not depend on time. The slope decreases for larger values of τ im . As described in detail in section 2.5 the long-time asymptote t ≫ τ m of the density follows a Gaussian with an effective advection speed and an effective diffusivity.

Mean squared displacement
In this section we analyse the MSD of the total density. In appendices E.1.1 and E.1.2 we analyse the MSD of the immobile and mobile density, respectively. We now investigate how advection changes the MSD as compared to the results presented in [54]. Table C1 shows a series of the total, mobile and immobile MSDs for initially mobile and initially immobile tracers for t ≪ τ m , τ im . In all cases, the leading order term does not depend on v. Therefore, for short times the MSDs are equivalent to the MSD without advection. In figure 6(a) the MSD for initially mobile tracers is shown for long immobilisations, τ m ≪ τ im , where the solid black line corresponds to the total density's MSD. In panel (b) we show the MSDs for short immobilisations, τ im ≪ τ m . For comparison, we show the MSD for the case without advection in grey. It can be seen that the corresponding MSDs with and without advection overlap for t ≪ τ v . After a linear  short-time behaviour for t ≪ τ ⋆ the MSD crosses over to a cubic scaling for τ ⋆ ≪ t ≪ τ m , τ im . Our goal is to quantify the cubic scaling and determine the value of τ ⋆ . From a series expansion of the exact expression for the MSD (exact expressions for the moments are given in appendix B), we obtain the asymptotic MSD In figure E2 we compare equation (27) to the full expression of the MSD and find very nice agreement. A cubic scaling of the MSD (27) emerges at intermediate times when the cubic term dominates over the linear term in equation (27). This corresponds to the relation Indeed, the cubic scaling of the total MSD is shown in figure 6(a) in the domain τ ⋆ ≪ t ≪ τ m . Notably, it occurs also for τ im ≪ τ m , as shown in figure 6(b). This cubic scaling is new compared to the case without advection. It lies within the domain where we found an advection dominated regime τ v ≪ t ≪ τ m in the previous section, in which the density (21) consists of a uniform and a Gaussian distribution with time-dependent weights. In figure E1 we show the MSD and choose such parameters that emphasise that the intermediate cubic scaling can be more than a bare crossover but corresponds to a distinct anomalous regime. We now show that the anomalous diffusion arises from this distribution. First we consider the uniform distribution ranging from zero to vt. It has the first moment ⟨x(t)⟩ im = vt 2 and the second moment 12 . This is the dominating term of the immobile density's MSD for τ v ≪ t ≪ τ m , τ im , as shown in figures 6(a) and (b), where the quadratic scaling can be observed. Second we recall the first and second moments of the mobile Gaussian distribution ⟨x(t)⟩ m = vt and ⟨x(t) 2 ⟩ = v 2 t 2 + 2Dt. Now we consider the MSD of the total density by combining the moments with the normalisations t/τ m and 1 − t/τ m for the uniform and Gaussian distribution, respectively. This leads to the same asymptotic MSD (27), as obtained by the series expansion, as the following calculation shows: The asymptotic expression is identical to what we found from the exact expressions (27). This shows indeed that the uniform distribution and the Gaussian distribution can indeed explain the cubic scaling. The cubic scaling of the MSD emerges for τ ⋆ ≪ t ≪ τ m , τ im . This means that for long immobilisations τ m ≪ τ im the advection needs to be sufficiently large such that τ v = 2D/v 2 ≪ τ m , as shown in figure 6(a). The cubic scaling emerges for short immobilisations τ im ≪ τ m , as well. In that case the Péclet number needs to satisfy Pe ≫ 3τ 2 m /τ 2 im . In appendix I we discuss the parameter regimes for which cubic scaling emerges in detail. Now we choose long immobilisations, τ im ≫ τ m , and consider the intermediate immobilisation dominated regime τ m ≪ t ≪ τ im . In the absence of advection we found a plateau of the total MSD shown as the grey solid line in figure 6(a). Compared to [54] we here choose a high Péclet number, v 2 τ m /2D ≫ 1, and keep the time scale separation τ m ≪ τ im here. The plateau still exists at τ m ≪ t ≪ τ im , although at a higher value as shown by the black solid line in figure 6(a). The existence of the plateau comes as no surprise, because the physical mechanism of the plateau remains unchanged. All tracers are initially mobile and have immobilised for τ m ≪ t ≪ τ im , as shown in figure 4 for t = 10. When all tracers are immobile the density does not change and hence the MSD remains constant. Analytically, this can be seen most easily from the expressions of the moments in Laplace space that we obtain from expression (A.5). For τ m ≪ t ≪ τ im , we have ϕ(s) ∼ 1/τ m , which is a constant. Since the Laplace inverse L −1 [1/s] = 1 we obtain a constant MSD in this time-domain τ m ≪ t ≪ τ im , where we observe the asymmetric Laplace distribution (23). We emphasise that the plateau exists regardless of the values of v and D. At long times t ≫ τ m , τ im , τ v the MSD grows linearly with the effective diffusion coefficient (τm+τ im ) 3 , as described in section 2.5.

Initially immobile tracers
In this section we consider initially immobile tracers. Experimentally, this may correspond to the situation when tracers are released into a microfluidic setup and (part of them) allowed to bind to the sensor receptors. Subsequently, the mobile tracers are flushed out, and then the recording is started. In section 4.1 we report the density for long immobilisations, τ m ≪ τ im , and high Péclet numbers. In section 4.2 we repeat the same steps for short immobilisations, τ im ≪ τ m . Section 4.3 is concerned with the MSD, and in the

Density for long immobilisations
We assume long immobilisations, τ m ≪ τ im , corresponding to the case studied in [54]. We choose a high Péclet number to emphasise the effect of advection, in contrast to the v = 0 case in [54]. At short times t ≪ τ v the propagator for mobile tracers is the same Gaussian as for the case v = 0. Therefore, we obtain the same short-time densities as in the v = 0 case, i.e. a δ-peak at the origin and an additional non-Gaussian distribution. At short to intermediate times, t ≪ τ m , τ im , at which t can be shorter or longer than τ v , most initially immobile tracers are concentrated at the origin and only gradually mobilise. This gives rise to the immobile density n im (x, t) ∼ (1 − t/τ im )δ(x). As described in detail in [54], immobile tracers that were initially mobile follow the same density as mobile tracers that were initially immobile, equation (19), up to a factor τ m /τ im . Therefore, we arrive at the total density for short to intermediate times for t ≪ τ m , τ im . The asymptote of the second summand corresponding to the mobile density in expression (31) is shown in figure 7 as the blue dashed line, which nicely matches the simulations and the Laplace inversions for t = 10 −3 and t = 10 −1 . We note that n tot (x, t) is always normalised, ∞ −∞ n tot (x, t)dx = 1, by construction. The same arguments as presented for initially mobile tracers explain the uniform density that appears at intermediate time scales for τ v ≪ t ≪ τ m , τ im , corresponding to the second panel in figure 7. In the immobilisation dominated intermediate time domain τ m ≪ t ≪ τ im the Laplace distribution with additional δ-peak emerges with the same scale parameter as for initially mobile tracers (23). The prefactor of the asymmetric Laplace distribution is now t/τ im , and this asymmetric Laplace distribution is shown in figure 7 as the grey dashed line. The long-time limit does not depend on the initial conditions and follows the same density as the initially mobile tracers, equation (17), as shown by the red dashed line in figure 7 at t = 10 4 .

Density for short immobilisations
We now consider short immobilisations, τ m ≪ τ im . At short to intermediate time scales t ≪ τ im ≪ τ m , at which t can be shorter or longer than τ v , the same expression (31) holds as for long immobilisations. This can be seen in figure 8, where the asymptote (31) is shown as the blue dashed line. We find excellent agreement between the histogram based on the simulations and the density (31) for t ≪ τ im , τ m , see the first and second panels in figure 8 for t = 2 × 10 −6 and t = 3 × 10 −4 , respectively. In contrast to the case v = 0 in [54], a new regime τ im ≪ t ≪ τ m emerges, that we call advection induced subdiffusion. The total density follows an exponential distribution, as can be seen from the linear shape of the density in the semi-log plot for t = 10 −2 . To explain this shape we solve the model equations (3) for the advection induced subdiffusion regime times t ≪ τ m , where t can be longer or shorter than τ v . This produces the total density where the integral is identical to expression (25) for immobile tracers that were initially mobile. We obtain the expression for v 2 τ im ≫ D and t ≪ τ m . This comes as no surprise, as in both cases we have the same Gaussian propagator, in which the tracers have varying immobile durations. In the case of initially mobile tracers the immobile duration stems from an immobilisation at t > 0. In the present case of initially immobile tracers the immobile duration arises from the slow release at the origin at t = 0. In both cases the immobile duration is drawn from an exponential distribution with mean τ im . The first term in (34) corresponds to initially immobile tracers that have not mobilised up to time t. The second term accounts for the slow release with rate τ −1 im exp(−t/τ im ) and motion in the mobile zone with advection only. The exponential distribution (34) is shown in figure 8 at t = 2 × 10 −2 as the grey dashed line, and we find good agreement with the Laplace inversion of n tot (x, s) and simulations.

Mean squared displacement
We now analyse the MSD of the total density. In appendices E.2.1 and E.2.2 we analyse the MSD of the immobile and mobile density, respectively. We show the MSDs for long and short immobilisations in panels (a) and (b) of figure 9, respectively. From a series expansion of the MSD at t = 0 we obtain the asymptotic MSD at short to intermediate times, in which t can be smaller or larger than τ v . We compare the asymptote (35) to the full expression of the MSD in figure E2(b) and find very nice agreement. At short times t ≪ τ v , the quadratic term dominates and we observe the same MSD as in the case without advection. This can be seen in figures 9(a) and (b). At intermediate times τ v ≪ t ≪ τ m , τ im the cubic term in the asymptotic MSD (35) dominates. As shown in figure 9 the cubic scaling emerges for short and long immobilisations. Now we go to the case of short immobilisation, τ im ≪ τ m , for which advection induced subdiffusion emerges for τ im ≪ t ≪ τ m . As described in section 4.2, the total density follows an exponential distribution (34) and a δ-peak at the origin for τ v ≪ t ≪ τ m . The MSD of that distribution is given by for τ v ≪ t ≪ τ m and τ im ≪ τ m , which is shown in figure 9(b) as the blue line. From expression (36) we recover the intermediate time asymptote implying the same cubic scaling for τ v ≪ t ≪ τ im as we found from the series expansion of the full MSD (35) for the advection dominated regime. This can be seen in the MSD in figure 9(b). The cubic scaling of the MSD emerges in the domain τ v ≪ t ≪ τ m , τ im , which limits the parameters to τ v ≪ τ m , τ im . In terms of the Péclet number this implies Pe ≫ 1 for long immobilisations τ m ≪ τ im and Pe ≫ τ m /τ im for short immobilisations τ im ≪ τ m . A detailed discussion of the parameter regimes and the coexistence of the plateau regime is presented in appendix J.
In the advection induced subdiffusion domain the MSD (36) reaches the plateau value These two anomalous scaling regimes shown in figure 9(b), namely, the cubic scaling and the plateau behaviour can be explained as follows. For τ v ≪ t ≪ τ im the MSD grows due to the slow release and fast advection, where the spread due to diffusion is negligible. When all tracers mobilised at t ≫ τ im , this spread due to advection vanishes, and the distribution moves along the direction of advection without changing the Table 1. Main results of anomalous scaling of the MSD. In the second column the short-time asymptotes of the total density's MSD is shown for initially mobile and initially immobile tracers. The cubic term dominates for the time regime given in the third column, with τv = 2D/v 2 and τ⋆ = √ 3τvτm. In the fourth column the time regimes of the plateaus are given. For initially immobile tracers this advection induced subdiffusion occurs for large Péclet numbers only. tracers initially t ≪ τm, τ im cubic regime plateau regime shape significantly. This explains the plateau in the MSD for τ im ≪ t ≪ τ m . We highlight that this advection induced subdiffusion is a new behaviour compared to the case without advection. For comparison, we show the MSD for v = 0 in figure 9(b) as the grey solid line. It crosses over from the short-time scaling ∼ Dt 2 /τ im to the long-time asymptote ∼ 2D eff t without any intermediate regime.

Conclusion
We analysed the densities along with the first and second moments of the MIM with exponential Poissonian switching between the mobile and immobile states in the presence of a drift velocity v. The whole dynamic is characterised by the mean mobile duration τ m , the mean immobile duration τ im and the time scale τ v = 2D/v 2 , which is related to the Péclet number Pe = τ m /τ v . For t ≪ τ v advection plays a negligible role and the process is diffusion dominated, yielding the same results as in [54], where the diffusion regimes were analysed for the advection-free case. In order to highlight the role of advection we choose τ v ≪ τ im , τ m , for which an advection dominated regime emerges for τ v ≪ t ≪ τ m , τ im . Relatively high Péclet numbers can be achieved in microfluidic setups and in certain geophysical systems. The first moment is proportional to the second moment of the advection-free model, as shown by the second Einstein relation. The second moment of the advection-free model has been discussed in detail in [54]. Therefore, we here concentrated on the discussion of the MSD. In general, for any fraction of initially mobile tracers and an arbitrary fraction τ im /τ m , we found the same long-time behaviour for t ≫ τ m , τ im . The total density follows a Gaussian with an effective advection speed v eff = vτ m /(τ m + τ im ) and an effective diffusivity D eff = D τm τm+τ im + v 2 τ 2 m τ 2 im (τm+τ im ) 3 . Compared to the advection v in the free phase, v eff is always smaller. The effective diffusivity D eff is always larger than the effective diffusivity in the advection-free case due to the velocity-dependent term. Specifically, for sufficiently high Péclet numbers the effective diffusivity can significantly exceed the diffusivity D in the mobile domain. While D eff was reported before [7,13], we here provided a physical explanation for the additional dispersion due to the variance of durations the tracers spent in the mobile state.
We analysed two specific initial conditions with fully mobile and fully immobile tracers in detail. For initially mobile tracers the advection dominated regime contains two parts. In the first part, 3τ v τ m , the MSD grows linearly in time. In the second part of the advection-dominated regime, τ ⋆ ≪ t ≪ τ m , τ im , the MSD grows cubically. This regime is valid for any ratio τ im /τ m . This can be seen in table 1, where we summarise the main results of the anomalous scaling of the MSD.
For long immobilisations τ im ≫ τ m for initially mobile tracers, and in the advection dominated regime, the density consists of a Gaussian shape of mobile tracers plus a spatially uniform distribution of immobilised tracers. At longer times, τ m ≪ t ≪ τ im , we recovered a plateau in the MSD of initially mobile tracers reported in [54], regardless of the presence of advection. In this immobilisation dominated regime the distribution follows an asymmetric Laplace distribution, which decays rapidly in the direction opposite of the advection velocity and falls of slowly in the direction of the advection velocity. In the opposite case of short immobilisations, τ im ≪ τ m , the total density follows a decreasing Gaussian shape with an increasing exponential tail in the direction opposite of the advection velocity in the whole advection dominated regime τ v ≪ t ≪ τ im . This exponential tail is similar to what was found for CTRWs with advection [46].
For initially immobile tracers the advection-dominated regime emerges, as well. This can be seen in table 1, where the series expansion of the MSD is shown. Here, the cubic scaling of the MSD appears in the whole regime τ v ≪ t ≪ τ m , τ im for any ratio of τ m /τ im and a sufficiently high Péclet number. This is in contrast to the advection-free case reported in [54], where the MSD is close to the Brownian case for short immobilisations, τ im ≪ τ m . For such short immobilisations the density follows a growing and drifting exponential distribution with an additional peak of immobile tracers at the origin. For later times, we found an advection induced subdiffusion regime τ im ≪ t ≪ τ m for high Péclet numbers and short immobilisations. This regime has not been reported previously. Here, the density consists of an exponential distribution with fixed scale parameters, where the mean is moving at a constant speed. For long immobilisations τ m ≪ τ im , the density follows a spatially uniform distribution for τ v ≪ t ≪ τ m . For completeness, we consider an equilibrium fraction of initially mobile tracers in, which naturally occur in experiments, in appendix K, where a ballistic scaling of the MSD emerges at intermediate times.
Now we put our work into some context. A CTRW with an exponential sojourn time distribution with mean τ and a Gaussian displacement density with non-zero mean µ and variance σ 2 was considered in [46]. This model may appear similar to our model with exponentially distributed residence times in the immobile state and advection-diffusion in the mobile state. However, the MSD ⟨[x(t) − ⟨x(t)⟩] 2 ⟩ = σ 2 t/τ + µ 2 t/τ is linear at all times for the CTRW, as shown in appendix L.1. This contrasts the anomalous scaling at intermediate times of the MSD found here. Another formulation of CTRW is the two-state CTRW, in which two transition densities are considered from which steps are drawn in an alternating way [51]. In [51] the long-time asymptotic diffusion coefficient is obtained, which matches our result for a specific choice of parameters. We also mention the case in [52] in which a CTRW is analysed, whose waiting time distribution function is the weighted sum of two exponentials. If one of the exponential distributions has a weight close to unity and a mean that is significantly shorter than τ m , this model generates the same MSD as the MSD of the total density from the MIM. We stress that the MIM provides additional information in the form of separate mobile and immobile densities. Thus the MIM is not simply the long-space-time limit of the model in [52]. A detailed analysis of the similarities and differences between the MIM and the processes in [51,52] deserves a separate study.
In contrast to the advection-free case of MIM, where intermittent anomalous diffusion solely arises for mobile tracers for long immobilisations, the anomalous transport behaviours reported here occur for both long and short immobilisations. A condition for this anomalous behaviour is a high Péclet number, which occurs for example in experiments with biomolecules in biosensors [38]. Short immobilisations may occur unintentionally due to unwanted binding to the surface of a flow cell or other experimental boundaries. The resulting transport will be anomalous at intermediate time scales, even if the tracers are only subject to Brownian motion with drift in the mobile state. While the anomalous diffusion and non-Gaussian densities occur at intermediate time scales only, the measurement time is finite in experiments, and the effects may therefore erroneously appear to be an asymptotic phenomenon. In conclusion, we found a variety of anomalous diffusion regimes and non-Gaussian displacement distributions at relevant intermediate time scales in a simple MIM. We emphasise that the model's simplicity is attributed to its dependency on merely two parameters in its dimensionless representation. We finally note that persistent-intermittent anomalous scaling behaviours of the MSD as reported here for seemingly simple Poissonian mobile-immobile dynamics may be relevant for numerous experimental settings. Such scenarios should therefore be included in contemporary data analysis methods by classical observables and machine learning approaches [67][68][69].

Data availability statement
No new data were created or analysed in this study.

Acknowledgments
We acknowledge funding from the German Science Foundation (DFG, Grant No. ME 1535/12-1). AVC acknowledges the support of the Polish National Agency for Academic Exchange (NAWA).

Appendix A. Solving the MIM model in Fourier-Laplace space
We apply the Fourier-Laplace transform to equations (3) and obtain the expressions as well as Fourier inversion directly produces expressions (4) to (6). From the densities (A.2) to (A.4) we obtain the pth moment (p ∈ N) in Laplace space. In order to obtain the mobile and immobile moments, we first calculate the non-normalised moments with j ∈ {m, im, tot}. We then normalise the moment (A.6) in the time-domain with the fractions of mobile and immobile densities,

Appendix B. Expressions of the moments
In this section we present the exact expressions of the first and second moments, which are obtained from expression (A.6).

B.1. First moments
For mobile initial conditions we find and For initially immobile tracers we find and

B.2. Second moments
Consider initially mobile tracers (f 0 m = 1 and f 0 im = 0). Then we obtain and For initially immobile tracers the moments read (B.12)

Appendix C. Series expansion of moments
The series expansion of the MSDs is shown in table C1.

Appendix D. Details on subordination
The probability density function P(τ, t) of τ (t) =´t 0 i(t ′ )dt ′ is calculated in [23] and reads Here, P u q,r (τ, t) and P b q,r (τ, t) denote the probability of initially mobile and immobile tracers, respectively, to spend in total the duration τ mobile in q mobile periods. The tracer remains immobile for the total duration t − τ in r immobile periods. The expressions for P u q,r (τ, t) and P b q,r (τ, t) are given by and In geological experiments, typically the mobile density is measured [18]. Similarly to expression (11), we obtain the mobile density with the probability P m (t ′ , t) to be mobile for a total duration t ′ at time t conditioned to be in the mobile phase at time t. From expression (D.1) we simply need to remove the terms with an additional immobilisation to obtain The immobile density n im (x, y, t) and P im (s, t) can be defined analogously to expressions (D.6) and (D.7) with The subordinator P(s, t) can also be obtained in another way, as we show now. Since i(t) follows a telegraph process, it is easy to see that the probability density functions P m (τ, t) ≡ P(τ, t|i(t) = 1) and P im (τ, t) ≡ P(τ, t|i(t) = 0) for the operational time obey the following set of equations where P m (τ, t) and P im (τ, t) denote the subordinator P(τ, t) conditioned on the tracer being mobile or immobile at time t, respectively, with P(τ, t) = P m (τ, t)

E.1. Mobile initial conditions
We assume initially mobile tracers and obtain asymptotic expressions in sections appendices E.1.1 and E.1.2.

E.1.1. MSD of the immobile density
Now we turn to the MSD of the immobile density. A series expansion of the exact MSD produces Expression (E.1) implies a ballistic scaling of the immobile MSD for 6τ v ≪ t ≪ τ m , τ im . This matches the MSD of the uniform distribution (20) of n im (x, t) shown in figure 4 where the right border moves at a constant speed. The immobile MSD with the ballistic scaling is shown in figure 6(a) as the dotted line. Panel (c) depicts the anomalous diffusion exponent α, which is close to two for 6τ v ≪ t ≪ τ m . Figure 6(a) shows the MSD of the mobile density as a dashed black line. For almost all times t ≪ τ m it coincides with 2 Dt corresponding to the case without advection and Brownian diffusion. The reason for this is that at t < τ m ≪ τ im almost all mobile tracers have never immobilised, as shown in figure 4.

E.1.2. MSD of the mobile density
The MSD of the mobile density has the series expansion The short-time behaviour 2 Dt dominates up to (2Dτ m τ im /v 2 ) 1/3 , as shown in figures 6(a) and (b). The time scale (2Dτ m τ im /v 2 ) 1/3 is very close to τ m in panel (a) and close to τ im in panel (b), therefore the t 4 scaling is not observed in figure 6. In the appendix in figure E1(a), we choose a higher Péclet number, for which the t 4 scaling is distinct. After τ m the MSD has a peak and decreases afterwards for some time. This can be explained similarly to the advection free case in [54]. Up to τ m , most of the mobile density consists of mobile tracers that have never immobilised. This can be seen in the histogram for t = 0.1 in figure 4, where most of the mobile density follows a Gaussian which is coloured black corresponding to zero immobilisation events N im = 0. When the mobile density mostly consists of tracers that were immobile once the MSD decreases because the leading Gaussian peak is missing.

E.2. Immobile initial conditions
We assume initially immobile tracers and obtain asymptotic expressions in sections appendices E.2.1 and E.2.2.
where t can be shorter or longer than τ v . In figure E2(b) we verify this asymptote and it shows good agreement with the full analytical solution. For a high Péclet number Pe = v 2 τm 2D ≫ 1 the quartic term in the asymptotic MSD (E.3) dominates for τ v ≪ t ≪ τ m , τ im , as shown in figures 9(a) and (b).

E.2.2. MSD of the mobile density
The MSD of the mobile density with immobile initial conditions is the same as the MSD (E.1) of the immobile density with mobile initial conditions [54]. In the long-time limit t ≫ τ v , τ m , τ im , the MSD is linear with effective diffusion coefficient (15). which corresponds in time-domain to the integral The integral (F.2) has the physical interpretation of mobile tracers following a Gaussian and immobilising with the constant rate 1/τ m .

Appendix G. Density with short immobilisations
For t ≪ τ m , the fraction t/τ m of mobile tracers is trapped for a short period τ drawn from γ(τ ) = exp(−τ /τ im )/τ im . This is shown as the red area in figure 5 denoting one immobilisation. Hence, these tracers were mobile for a total period of t − τ . We convolute this with the propagator for advection diffusion and obtain the expression with the complimentary error function erfc(x). Note that (G.1) is also valid for v 2 < 4D/τ im , i.e. also for the case without advection. The seemingly imaginary parts cancel out. Expression (G.1) can be rewritten as a(x, t) = Im(erf( ibt−|x| √ Dt ) exp(−ib|x|/2D)) − sin(b|x|/2D)/b, which is strictly positive with ib = v 2 − 4D/τ im , and it does not oscillate. The − sin(. . .) rather removes the oscillatory part of the Im(exp(−ib|x|/2D)) and the error function approaches unity for large x.

Appendix H. Dependence of D eff on the model parameters
We now analyse the dependence of the effective diffusion coefficient on the model parameters. In section 2.3 we show that only two dimensionless parameters characterise the model, namely the Péclet number and τ im /τ m . We use the decadic logarithms of these parameters as the abscissa and ordinate in figure H1. The logarithmic colour map encodes the ratio D eff /D for D eff > D. Outside that region D eff is smaller than D, as in the case without advection in [54]. For high Péclet numbers D eff takes on high values depending on τ im /τ m . Now we analyse this dependence in more detail. For D > 0 we now analyse the ratio D eff /D = c(Pe, τ m /τ im ), that depends on Pe and τ m /τ im . For c ⩾ 1 we find the expression which is the Péclet number as a function of τ im /τ m . It is shown for a range of c values in figure H1 as solid grey lines and matches the colour plot. We obtain the asymptotes  Depending on the parameters, the MSD of the total density with initially mobile tracers has a cubic scaling for τ⋆ ≪ τm, τ im and a plateau for τm ≪ t ≪ τ im . We quantify the existence of the cubic scaling with the maximum value of the anomalous diffusion exponent α for each point in the map. We show this αmax with a red colour coding. The black dotted line shows the border of the domain where the duration of the cubic scaling is larger than zero, i.e. t3 > 0. It has the asymptote ≃ Pe −1/2 for large Péclet numbers. The black dashed line shows the contour line, at which D eff = 100D. It has the same asymptotic scaling ≃ Pe −1/2 for large Péclet numbers and short immobilisations. For short immobilisations, i.e. negative ordinates, the shape of the read area follows the same asymptote. We quantify the existence of the plateau with the minimum value of the anomalous diffusion exponent for each point in the map. We show this α min with the blue colour coding. The top left corner corresponds to the parameters discussed in [54]. In the top right corner we have a coexistence of the cubic scaling and the plateau in the MSD. In the white area both αmax and α min are close to unity, i.e. the MSD is close to linear at all times.
τ im < τ m afterwards. Rearranging the condition τ ⋆ < τ m for the Péclet number reveals the lower bound 27/16 < Pe independent of τ im /τ m . Indeed, the black dotted line for t 3 = 0 is vertical in the map (I1) for τ im ≫ τ m , indicating no dependence on τ im /τ m in that regime. So far, we looked at the necessary condition that the lower bound for the interval of the cubic regime is lower than the upper bound. Another approach is to analyse the instantaneous anomalous diffusion exponent α(t). It can be obtained by calculating the slope in a double-logarithmic plot. This is shown in figure E1(c), where α remains close to unity for t < τ ⋆ and reaches the value three for τ ⋆ ≪ t ≪ τ m , τ im . The alternative criterion for the cubic regime to exist is then simply that the maximum value of α max is close to three. With this definition of α(t), superdiffusive parameter regimes are shown with the red colour-coding ranging from one to three in figure I1. This means that the more intense the red in an area is, the closer the parameters are to the asymptotic limits in which the cubic regime emerges. We notice that the red area lies entirely in the region with t 3 > 0 and its border follows the same asymptote Pe −1/2 for τ m ≫ τ im . The cubic regime appears both for D eff > D and for D eff ≲ D for long immobilisations τ im ≫ τ m and for Pe ≫ 1, as shown in the top right corner of the map in figure I1. For reference, a contour plot of D eff /D is shown in figure H1. The reason for the appearance of the cubic regime regardless of the ratio D eff /D is the existence of the plateau for τ im ≫ τ m , which we explain now in detail. We know that regardless of the ratio τ im /τ m the MSD of initially mobile tracers has the asymptote 2 Dt for short times t ≪ τ m , τ im , τ v . In the long-time limit it has the asymptote 2D eff t. The MSD of the total density is a strictly monotonic function. This means that if D eff < D, any growth faster than linear must be compensated by a regime with a sublinear growth. This slower than linear growth arises for long immobilisations only and constitutes the plateau for a sufficiently large ratio τ im /τ m .

I.3. Short immobilisations
Finally, we analyse the case of short immobilisations, τ im ≪ τ m . The same series expansion of the total MSD holds as in the case with long immobilisations (27). The difference between the upper and lower bound of the cubic scaling is therefore t 3 = τ im − τ ⋆ . Notably, the duration of the cubic regime is shorter as compared to the case with long immobilisations in figure 6(a), because τ ⋆ = 6Dτ m /v 2 is close to the lower bound τ im for τ m ≫ τ im . We emphasise that the cubic regime appears for short immobilisations, whereas for the case without advection the MSD is close to normal, as shown by the grey line in figure 6(b). For τ im ≪ τ m , i.e. for negative ordinates in figure H1, the condition τ ⋆ ≪ τ im simplifies to D v 2 τm = (2Pe) −1/2 ≪ 2/3 τ im τm . This means that for each ratio τ im /τ m there exists a lower bound for the Péclet number for the cubic regime to exist. Indeed, the black dotted line on the map in figure I1 for t 3 > 0 has the asymptotic scaling Pe −1/2 for τ m ≫ τ im . The cubic regime appears for t 3 > 0 only, i.e. only in the top right part from the black dashed line in figure I1. In figure I1 the red colour denotes the maximum anomalous diffusion coefficient α max . The condition D eff > 100D has the same asymptotic regime Pe −1/2 for τ m ≫ τ im , as shown by the black dashed line in figure I1. For τ m ≫ τ im the dark red region with α max close to three is almost entirely contained within this region. This means that for τ m ≫ τ im the cubic regime always coexists with D eff ≫ D. In contrast, we found that the cubic regime appears regardless of the ratio D eff /D for long immobilisations.
In addition to the maximum value of the instantaneous anomalous diffusion exponent, we consider the minimum value α min in figure I1, in which values close to zero correspond to very slow intermediate growth, i.e. a plateau. This allows us to analyse how the plateau at τ m ≪ t ≪ τ im is influenced by the presence of advection and the cubic regime. The coexistence of the cubic regime and the plateau can be seen in the MSD in figure 6(a). The top left blue corner in figure I1 corresponds to the case considered in [54] with long immobilisations and (almost) no advection. The blue subdiffusive area does not significantly change for growing Péclet numbers and it overlaps with the red superdiffusive region in the top right corner of the map in figure I1. In the lower left part of figure I1(a) white area dominates, in which diffusion is close to normal, with both α min and α max close to unity. This can be seen from the solid grey line in figure 6(b). Figure J1. Map with sub-and superdiffusive regimes, same as figure I1 for initially immobile tracers. The y-axis is chosen as log 10 (τ im /τm) and is an indicator for the strength of the effect of the immobilisations on the MSD. The x-axis is the logarithm of the Péclet number Pe = v 2 τm/2D. Both parameters span the full phase space, as described in section 2.3. Depending on the parameters, the MSD of the total density with initially mobile tracers has a cubic scaling for τv = v 2 /2D ≪ τm, τ im and advection induced subdiffusion for τ im ≪ t ≪ τm. We quantify the existence of the cubic scaling with the maximum value of the anomalous diffusion exponent α for each point in the map. We show this αmax with a red colour coding. The black dashed line shows the scaling ≃ Pe −1 for large Péclet numbers. We quantify the existence of the advection induced subdiffusion by the minimum value of the anomalous diffusion exponent for each point in the map. We show this α min with the blue colour coding. The top left corner corresponds to the parameters discussed in [54]. In the bottom right corner we have a coexistence of the cubic scaling and the plateau in the MSD. In the white area αmax is close to two and α min is close to unity. The white contour lines for αmax = 2.5 and α min = 0.5 serve as a guide to the eye. Figure K1. Same as figure 6 for equilibrium initial conditions f 0 m = f eq m = τm/(τm + τ im ). In contrast to the case with v = 0 in [54], the MSD of the total density (solid lines in (a) and (b)) are no longer linear at all times. A quadratic scaling emerges at τv/f eq im ≪ t ≪ τm, τ im . The MSD of the mobile density is identical to the MSD of the total density with initially mobile tracers in figure 6.

Appendix L. Comparison to CTRW
In this section we compare our results to three versions of CTRWs. For didactic purposes, we start in appendix L.1 with the classical CTRW with exponentially distributed waiting times. In appendix L.2 we compare our model to the two state CTRW [51]. Finally, in appendix L.3 we compare our model to a CTRW in which the waiting times are drawn from a waiting time distribution containing two exponential distributions [52].

L.1. Classical CTRW
We here compare the results obtained from the MIM (3) to a CTRW with advection and exponentially distributed sojourn times. The relation between MIM and CTRW has been discussed in detail, for example, in [13,48,70]. Let us define the sojourn time density ψ(t) = exp(−t/τ )/τ and the Gaussian displacement density λ(x) with mean µ and variance σ 2 . This corresponds to the model considered in [46]. In [71] the solution of the density function p(x, t) in this CTRW framework is given as for µ = 0. We incorporate the non-zero mean µ and obtain the solution Let us turn to the moments of p(x, t) (L.2). The fist moment is given by (L.3) Figure L1. Comparison of the MSD from MIM and CTRW with double exponential waiting time [52]. CTRW data originate from simulations and the full expressions from B are used to evaluate the MSDs of MIM. Figure L2. Mean number of steps and variance of step numbers for CTRW with double exponential waiting time [52] and MIM. CTRW data originate from simulations and the full expressions from appendix B are used to evaluate the MSDs of MIM. It is seen that both models shows perfect agreement in the drift-free case, for the specific choice of parameters indicated after equation (L.9), whereas in the presence of a drift, the models converge only in the long time limit.
of steps E(n) (calculated by setting D = 0), where the dots for CTRW and the line for MIM overlap at all times. Therefore, the MSDs in the advection-free case of CTRW and MIM in figure L1 are equivalent. In contrast, in the case of advection, the variance of step numbers var(n) couples to the MSD, as can be seen in expression (L.10). The variance of step numbers differs for t < τ D . We explain this as follows. If a tracer is in the mobile state of MIM, it will have a fixed number of steps n = t/∆t at time t. The stochasticity of step numbers arises solely from immobilisations. By setting D = 0 we obtain a cubic short-time growth of the variance from the short-time asymptote (27). In contrast, in the CTRW case, the number of steps is random even in the case of a simple CTRW with only one exponential waiting time distribution, as described in appendix L.1. It can be seen in the variance of CTRW (L.5) for σ = 0 and µ = 1, where the variance of step numbers grows linearly. This is due to the fact that after each step a random waiting time is drawn from the exponential distribution. Therefore, the variance of step numbers is higher for small times and hence the MSDs of MIM and CTRW do not overlap for t < τ D , as shown in figure L1. We point out that the disagreement between the MSD of MIM and CTRW with advection lasts until ≈ 3 × 10 −1 , which is notably larger than the chosen τ D = 2 × 10 −2 .