Tracer Diffusion on a Crowded Random Manhattan Lattice

We study by extensive numerical simulations the dynamics of a hard-core tracer particle (TP) in presence of two competing types of disorder - frozen convection flows on a square random Manhattan lattice and a crowded dynamical environment formed by a lattice gas of mobile hard-core particles. The latter perform lattice random walks, constrained by a single-occupancy condition of each lattice site, and are either insensitive to random flows (model A) or choose the jump directions as dictated by the local directionality of bonds of the random Manhattan lattice (model B). We focus on the TP disorder-averaged mean-squared displacement, (which shows a super-diffusive behaviour $\sim t^{4/3}$, $t$ being time, in all the cases studied here), on higher moments of the TP displacement, and on the probability distribution of the TP position $X$ along the $x$-axis. Our analysis evidences that in absence of the lattice gas particles the latter has a Gaussian central part $\sim \exp(- u^2)$, where $u = X/t^{2/3}$, and exhibits slower-than-Gaussian tails $\sim \exp(-|u|^{4/3})$ for sufficiently large $t$ and $u$. Numerical data convincingly demonstrate that in presence of a crowded environment the central Gaussian part and non-Gaussian tails of the distribution persist for both models.

Quenched (frozen) spatial disorder which entails a temporal trapping of a tracer particle (TP) at some positions, often produces an anomalous sub-diffusive behaviour, especially in low-dimensional systems. Here, the TP trajectories are spatially more confined than the trajectories of a standard Brownian motion. As a consequence, the disorder-averaged mean-squared displacement (DA MSD) behaves as < X 2 (t) >∼ t γ , with t being time and γ -the dynamical exponent which is less than unity. Here and henceforth, the bar denotes averaging over thermal histories while the angle brackets stand for averaging over disorder. Striking examples of such a dynamical behaviour are provided by, e.g., the so-called Sinai diffusion in one-dimensional systems [14] (see also Refs. [3][4][5][6]) in which the DA MSD grows as < X 2 (t) >∼ ln 4 t (i.e., formally, γ = 0), Sinai diffusion in presence of a constant external bias [15,16] or migration of excited states along a one-dimensional array of randomly placed donor centres [1,6]. In this latter example the dynamical exponent γ is non-universal and equals the mean density of donor centres times the characteristic length-scale of the distance-dependent (exponential) transfer rate. If this product is less than unity, a sub-diffusive motion takes place. Two other examples concern diffusion in the "impurity band" [17] and the so-called Random Trap model [18][19][20]. Here, as well, γ is non-universal and is less than unity in some region of the parameter space. In higher-dimensional systems, diffusion in presence of such a disorder typically becomes normal (see, however, Refs. [17,21,22]) and the disorder affects only the value of the diffusion coefficient. Diffusion is also normal in the asymptotic large-t limit in one-dimensional systems with a periodic disorder. Here, however, the value of the diffusion coefficient may exhibit strong sample-to-sample fluctuations and thus have non-trivial statistical properties, such that the averaged diffusion coefficient will not be representative of the actual behaviour (see, e.g., Ref. [23]). The large-t relaxation of the diffusion coefficient to its asymptotic value may shed some light on the kind of disorder one is dealing with [24].
Random frozen convection (velocity) flows most often produce a super-diffusion with γ > 1. To name just two such situations, we mention a model in which a TP is passively advected by quenched, layered, randomly-oriented flows (say, along the x-axis) and undergoes a normal diffusion in the direction perpendicular to them (i.e., along the y-axis), as well as its generalisation -a random Manhattan lattice (see Fig.  1), in which the orientation of convection flows randomly fluctuates both along the streets and avenues (i.e., along both x-and y-axes). The former model was introduced originally for the analysis of conductivity of inhomogeneous media in a strong magnetic field [25] and of the dynamics of solute in a stratified porous medium with flow parallel to the bedding [26]. In such a setting, usually referred to as the Matheron -de Marsily (MdM) model according to the names of authors of Ref. [26], the TP dynamics in the flow direction (along the x-axis) is characterised by a super-diffusive law of the form < X 2 (t) >∼ t 3/2 , i.e., γ = 3/2. Many interesting generalisations and more details on the available analytical and numerical results can be found in Refs. [27][28][29][30][31][32][33][34][35]. Diffusion of a single TP on a square random Manhattan lattice has been analysed in Refs. [27,28]. It was shown, by using simple analytical arguments and a numerical analysis, that in this case the DA MSD also exhibits a super-diffusive behaviour, but with a somewhat smaller dynamical exponent γ = 4/3, i.e., the DA MSD of the x-component of the TP position obeys < X 2 (t) >∼ t 4/3 . This model has been also widely studied in different contexts in mathematical literature (see, e.g., Ref. [36]). A generalisation of a random Manhattan lattice was invoked as an example of a plausible geometric disorder in a recent analysis of the localisation length exponent for plateau transition in quantum Hall effect [37]. This latter setting, however, is clearly more complicated than the MdM model with the layered flows and the theoretical progress here is rather limited; the behaviour beyond the temporal evolution of a DA MSD is still largely unknown.
Dynamical disorder emerges naturally when the TP's transition rates fluctuate randomly in time, as it happens, for instance, in physical processes underlying the socalled diffusing-diffusivity models [38][39][40][41][42][43][44] or the dynamic percolation [45][46][47]. Another pertinent case concerns the situations when the TP evolves in a dynamical environment of mobile steric obstacles -interacting crowders which impede its dynamics (see, e.g., Refs. [10][11][12]). A paradigmatic example of such a situation is provided by a TP diffusion in lattice gases of hard-core particles, which undergo the so-called simple exclusion process (see Ref. [13] for a recent review), i.e., perform lattice random walks subject to the constraint that each lattice site can be at most singly occupied. It is well-known that in such an environment the particles' dynamics is strongly correlated. These correlations are especially important and cause an essential departure from standard diffusive motion in two cases: a) in one-dimensional geometry -the socalled single-files, in which the particles cannot bypass each other and the initial order of particles is preserved at all times; and b) on ramified comb-like structures consisting of an infinitely long single-file backbone with infinitely long single-file side branches, which permit for some re-ordering of particles. In single-files, the TP meansquared displacement exhibits an anomalous sub-diffusive behaviour X 2 (t) ∼ t 1/2 . This striking result was first obtained analytically by Harris [48] (see Refs. [49,50] for a review), and holds also for all the cumulants of X(t) [51,52] and in case of multiple TPs [53][54][55]. On crowded comb-like structures, the TP mean-squared displacement exhibits a variety of sub-diffusive transients and, in some cases, an ultimate subdiffusive behaviour [56]. On higher-dimensional lattices, the TP dynamics becomes diffusive in the large-t limit with the effective diffusion coefficient being a non-trivial function of the density of crowders and other pertinent parameters [57][58][59][60][61][62][63]. This nontrivial behaviour of the diffusion coefficient is associated with the enhanced probability of backward jumps -in a crowded environment, for any particle it is more probable to return back to the site it just left vacant, than to keep on going farther away [57][58][59][60][61][62][63].
Meanwhile, a considerable knowledge is accumulated through case-by-case theoretical and numerical analyses of the TP dynamics in a variety of model systems with either quenched or dynamical disorder (see, e.g., Refs. [1][2][3][4][5][6][7][8][9][10][11][12][13] and references therein). On contrary, still little is known about the TP diffusion in situations in which several types of disorder are acting simultaneously. To the best of our knowledge, the only work addressing specifically this question is recent Ref. [50], which focused on the TP random motion in single-files of hard-core particles having a broad scale-free distribution of waiting times, e.g., due to a temporal trapping of particles. Using some subordination arguments and numerical analysis, it was shown that here a combined effect of the disorder in transition rates and of the dynamical environment leads to a severe slowing-down of the TP random motion. Namely, the DA MSD of the TP follows X 2 (t) ∼ ln 1/2 (t), i.e., exhibits an essentially slower growth with time than the one taking place in systems in which either type of disorder is present alone. In case when a characteristic mean waiting time exists, i.e., the distribution is not scalefree, but the second moment diverges, the DA MSD grows faster than logarithmically, X 2 (t) ∼ t γ with γ < 1/2, but still slower than the above mentioned Harris' law. This paper is devoted to a question of the TP dynamics in presence of two interspersed types of disorder, which act concurrently and compete with each other. We consider the TP random motion subject to quenched random convection flows, which prompt a super-diffusive behaviour of the TP, in a dynamical environment which is damping its random motion. More specifically, we study here by extensive numerical simulations the dynamics of a TP which evolves on a square random Manhattan lattice of frozen (i.e. not varying in time) convection flows in presence of a lattice gas (LG) of mobile hard-core particles. The latter are either insensitive to convection flows, performing standard random walks among the nearest-neighbouring sites of a lattice with the probability 1/4 to go in any direction (Model A), or follow the convection flows (similarly to the TP) by choosing randomly between the two directions prescribed by a local directionality of bonds of the random Manhattan lattice (Model B). In the latter case the backward jumps of any LG particle are completely suppressed. The backward jumps of the TP are forbidden in both models. For both models, the TP and the LG particles obey a simple exclusion constraint, which effectively correlates the TP random motion and the evolution of LG particles. We focus on such characteristics of the TP dynamics as its DA MSD, and generally, the moments of arbitrary order, the distribution of its position at time moment t averaged over disorder, as well as the time evolution of the kurtosis of this distribution. We also address a question of the sample-to-sample fluctuations and analyse the MSD of the TP and the probability distribution of its position for several fixed realisations of disorder.
The paper is outlined as follows: In Sec. 2 we define the model under study and introduce basic notations. In Sec. 3 we discuss dynamics of a single TP in absence of the LG particles, appropriately revisiting the arguments presented in Refs. [27,28]. We also present here results of numerical simulations for the DA MSD and for higher moments of the TP displacement, as well for the disorder-averaged probability distribution of the TP position along the x-axis. This sets an instructive framework for the analysis of the TP dynamics in presence of LG particle. We close Sec. 3 addressing the issue of sample-to-sample fluctuations and also examine the spectral properties of the TP trajectories, which reveal several interesting features. In Sec. 4 we consider the TP dynamics in presence of LG particles for both Model A and Model B. Finally, in Sec. 5 we conclude with a brief recapitulation of our results.

Model
Consider a two-dimensional random Manhattan lattice (see Fig. 1), i.e., an infinite in both directions square lattice with unit spacing, decorated with arrows in such a way that directionality of each of them is fixed along each street (East-West) and an avenue (North-South) for their entire length, but whose orientation varies randomly from a street (an avenue) to a street (an avenue).
Let an integer n, n ∈ (−∞, ∞), numerate the columns (avenues) of the lattice, and an integer m, m ∈ (−∞, ∞), -the rows (streets), respectively. Then, the pattern of arrows in a given frozen realisation of convection flows is specified by assigning to each lattice site (with integer coordinates (n, m)) a pair of quenched random, mutually uncorrelated "bias" variables η n and ζ m . We use a convention that η n = +1 if an arrow points to the North, and η n = −1, otherwise; and ζ m = +1 if an arrow points to the East, and ζ m = −1, otherwise. We focus solely on the case when there is no global bias; that being, η n and ζ m assume the values ±1 with equal probabilities, which implies that η n = ζ m = 0. Furthermore, we stipulate that there are no correlations between the directions of arrows at n and n ′ , and at m and m ′ , i.e., η n η n ′ = δ n,n ′ , where δ a,b is the Kronecker-delta, such that δ a,b = 1 for a = b, and equals zero otherwise.

A single tracer particle
At time moment t = 0 (t is a discrete time variable, t = 0, 1, 2, . . .), we introduce the TP at the origin of the lattice and let it move, at each tick of the clock, according to the following rules: -at each discrete time instant t, we toss a two-sided "coin" ξ t which can assume, (with equal probabilities = 1/2), the values +1 and −1.
-being at position R t = (X t , Y t ), (where X t and Y t are the projections of R t on the x-and y-axes), the TP is moved, after choosing the value of ξ t , to a new position where the vectorial increment δ t is defined as with e x and e x being the unit vectors in the x-and y-directions, respectively. The expression (2) can also be conveniently rewritten in form of two coupled, non-linear recursion relations for the integer-valued components X t and Y t : Therefore, once (with probability 1/2) ξ t = 1, the TP is moved onto the neighbouring site along the x-axis in the direction prescribed by ζ Yt , and does not change its position along the y-axis. Conversely, if ξ t = −1, the TP is moved on a unit distance along the y-axis in the direction prescribed by η Xt , and does not change its position along the x-axis. We recall that the ensuing motion of the TP as defined by the recursion relations (4) is super-diffusive, with the dynamical exponent γ = 4/3 [27,28]. We note parenthetically that it may apparently be possible to find an equivalent two-dimensional model in the continuum space and time limit, write down coupled Langevin equations for the time evolution of the components and, eventually, define the associated Fokker-Planck equation obeyed by the probability Π t (X, Y ) of finding the TP at position (X, Y ) at time moment t for a given realisation of disorder. We will address this question in our following work. Second, it was claimed in Refs. [27,28] that at a coarse-grained level the TP dynamics on a random Manhattan lattice becomes equivalent to a Brownian motion in continuum, in a divergenceless random velocity field with power-law decay of the velocity correlation function. We however remark that going to a continuum limit necessitates a generalisation of the model studied here; in our settings, the jumps against an arrow are not permitted which tacitly presumes that the force acting on the particle along a given bond is infinitely large. Therefore, one has to allow for the jumps against an arrow and let them occur with a smaller (but finite) probability, than the probability of the jumps along an arrow. This is tantamount to considering finite forces. We, however, do not expect any substantial change in the dynamics in the finite force case, as compared to our model.
The algorithm of our numerical simulations of the TP dynamics on a random Manhattan lattice follows the relations (4). We generate trajectories along the xabd y-axes of a given length t, for a given set of thermal variables {ξ t } and a given realisation of "bias" variables η n and ζ m . The obtained individual trajectories are stored and the characteristic properties of interest -the moments of the TP displacement and the distribution function of the TP position -are evaluated by averaging over different realisations of trajectories. Averaging is first performed over 10 4 trajectories generated for a fixed realisation of a random Manhattan lattice, and then the procedure is repeated for 2 × 10 5 realisations of disorder. Simulations are performed for lattices containing L×L sites with L = 2×10 6 . Care is taken that neither of the TP trajectories reaches the boundaries of the lattice within the observation time, such that the finite-size effects do not matter. For the lattice size used in our numerical modelling, this permits us to safely explore the TP dynamics for times up to t = 10 6 . Lastly, we also analyse the sample-to-sample fluctuations and, in particular, address a question of the TP dynamics in presence of a single fixed realisation of disorder. In this case, for a given random realisation of disorder we run 2 × 10 7 trajectories.

The TP dynamics on a crowded random Manhattan lattice
The TP dynamics on a random Manhattan lattice populated with N − 1 lattice gas particles is analysed numerically. Due to a significant number of the particles involved, we are only able to consider square lattices with the maximal linear extent L = 2×10 3 . This means that the maximal time t, until which the finite-size effects can be discarded, is of order of 4 × 10 3 . Moreover, due to computational limitations, we record only 50 TP trajectories for each given realisation of disorder, and average over 10 3 realisations of disorder. Such a statistical sample appears to be sufficiently large to probe the behaviour of the DA MSD of the TP, but does not permit us to make absolutely conclusive statements about the shape of the distribution function. Nonetheless, our numerical data rather convincingly demonstrate that the overall behaviour of the latter is very similar to the one observed for the TP dynamics in absence of the LG particles, in which case a more ample statistical analysis has been performed.
The simulations are performed as follows: We first place the TP at the origin of a lattice and then distribute N − 1 hard-core particles among the remaining sites by placing a LG particle at each lattice site, at random, with probability ρ = N/L 2 . The latter parameter defines the mean density of particles in the system; in our simulations, we study the TP dynamics for nine values of ρ, ρ = 0.1, 0.2, 0.3, . . . , 0.9.
After the particles are introduced into the system, they are let to move randomly subject to a single-occupancy constraint. We distinguish between two possible scenarios: In model A we suppose that all the LG particles are not sensitive to the frozen pattern of convection flows and perform symmetric random walks, subject to the constraint that there may be at most a single particle (i.e., either the TP or a LG particle) at each lattice site. On contrary, for the TP the choice of the jump direction is dictated by the arrows present at the site it occupies at time moment t. As described above, the TP chooses at random between the two arrows outgoing from the site it occupies. In this case, the TP (which exhibits a super-diffusive motion in absence of the LG particles) is not identical to the LG particles and moves in a quiescent "fluid" of hard-core particles which exerts some frictional force on it. Note that here the backward jumps are forbidden for the TP only.
More specifically, at each step we select at random a particle, which can be either a TP or a LG particle, and let it choose the jump direction: if the selected particle is a TP, it chooses at random between the two arrows. Conversely, a LG particle chooses at random one among four neighbouring sites with probability = 1/4. The jump of a TP or a LG particle is fulfilled, once the target site is empty at this time instant; otherwise, the particle remains at its position. The time t is increased by unity after repeating such a procedure N times, such that all N particles present in the system, on average, have a chance to change their positions.
We have already mentioned that in this model the dynamics of LG particles is rather non-trivial due to an enhanced probability of backward jumps; it means that a particle which jumps onto an empty target site will most likely return on the next time step to the site it just left vacant, then will keep on going away from it. Even in absence of the TP and random convection flows acting on it, this circumstance results in a non-trivial dependence of the self-diffusion coefficient D tp of any tagged particle on the overall density of the LG particles. This dependence is known only in an approximate form (see, e.g., Refs. [13] and [57][58][59][60][61][62][63]). The available exact results concern the leading, in the dense limit ρ ≃ 1, behaviour of the self-diffusion coefficient D tp ≃ (1 − ρ)/(4(π − 1)) [64] and of the mobility µ tp ≃ β(1 − ρ)/(4(π − 1)) [65] of a tagged particle subject to a vanishingly small external force, with β being the reciprocal temperature. The appearance of the Archimedes' irrational number "π" seems astonishing and points on a non-trivial behaviour.

Model B.
In this model, we suppose that all the particles in the system are identical. It means that both the TP and the LG particles move on the lattice subject to a single-occupancy constraint and obey the rules of the random Manhattan lattice, by following the jump directions prescribed by the arrows.
Note that in this model the backward jump probability is equal to zero for all the particles, both for the LG particles and the TP. As a consequence, we expect that here the environment in which the TP moves is a kind of a "turbulent" fluid, in which all the particles exhibit a super-diffusive motion. Hence, we may expect that the environment becomes perfectly stirred at sufficiently large times, such that the time t gets simply rescaled by the frequency (1 − ρ) of successful jump events, (which is not the case for Model A). We are going to verify if this is the case in what follows.
3. Dynamics of a single tracer particle

Disorder-averaged mean-squared displacement
In order to calculate the DA MSD of a single TP moving on a random Manhattan lattice in absence of the LG particles, we suitably revisit the arguments presented in Refs. [27,28]. The latter were based on an estimate of typical fluctuations of sums of quenched random variables η n and ζ m , and a plausible closure relation. Here, we pursue a bit different line of thought.
First, we "solve" the recursions in eqs. (4) to get, for t ≥ 1, with the initial condition X 0 = Y 0 = 0. Expressions (5) define the TP positions X t and Y t for any t, for fixed realisations of thermal noises ξ t and "biases" η n and ζ m . We concentrate on the x-component and write down formally its squared value: Let the bar denote averaging over ξ τ -s, which amounts to averaging over thermal histories, and the angle brackets -averaging over random variables η n and ζ m , i.e., averaging over quenched disorder. Consider the averaged first sum in the right-handside (rhs) of eq. (6). Noticing that ζ 2 Yτ ≡ 1, i.e., is not fluctuating, we realise that the averaged first sum is simply Hence, the contribution of the averaged first sum to the DA MSD of the TP along the x-axis is that of a standard, discrete-time random walk (with the diffusion coefficient D = 1/4) on a two-dimensional undecorated square lattice. Focus on the summand in the second term in the rhs of eq. (6) and write down formally its averaged value: Note that we are allowed to perform averaging over ξ τ ′ directly, due to the fact that both Y τ ′ and Y τ are statistically independent of a random variable Note that only the product ζ Yτ ζ Y τ ′ is dependent on random convection flows. Averaging this product over quenched disorder, we find that, in virtue of the definition in eq. (1), expression (8) takes the form i.e., it is an averaged over thermal noises product of the indicator functions of two As a consequence, the expression 9 is the joint probability P (Y τ |t = τ ′ ; Y τ |t = τ + 1; Y τ |t = τ ) of the events that the TP trajectory Y t , with Y t=0 = 0, a) paused at its (unspecified) position at t = τ and b) returned at time moment t = τ ′ to the position it occupied at t = τ . The probability P (Y τ |t = τ ′ ; Y τ |t = τ + 1; Y τ |t = τ ) decouples into the product of the probability that the trajectory Y t appeared at an unspecified position Y τ at time moment t = τ , which equals unity since averaging over ξ k with k ∈ (0, τ − 1) implies averaging over all possible Y τ ; the probability that Y t paused at t = τ , which equals 1/2; and the probability that Y t returned to Y τ = Y τ +1 within τ ′ − τ − 1 steps. Making a plausible assumption that the dynamics, at least in the asymptotic limit t → ∞, does not depend of the starting point, we thus find that the expression (9) reduces to where P τ ′ −τ −1 (Y = 0) is the probability that the y-component of the TP trajectory returns to Y = 0, (not necessarily for the first time), on the (τ ′ − τ − 1)-th step.
Here, P t (Y ) (P t (X)) is a marginal distribution obtained from the full probability distribution function P t (X, Y ) of finding the TP at site (X, Y ) at time moment t by summing the latter over all X (Y ), that is Summing up the presented above reasonings, we arrive at the following representation of the DA MSD:  (20) and (21)). A superdiffusive behaviour sets in from rather short times and the diffusive transient (see the first term in the rhs of eq. (12)) is not observed. The inset displays the rate of a convergence of the time-dependent dynamical exponent γt, eq. (18), to its asymptotic value 4/3. Panel (b). Reduced moments < |Xt| q > t 2q/3 as functions of time. The dashed lines (from top to bottom) correspond to m 4 = 1.038, m 3 = 0.687, m 1 = 0.590 and m 2 = 0.556 (see eq. (21)).
Further on, the probability P τ ′ −τ −1 (Y = 0) is evidently a decreasing function of the difference τ ′ − τ − 1. Very general arguments (see also the numerical results presented in Fig. 3, panel (a)), suggest that P τ ′ −τ −1 (Y = 0) decays as a power-law : in the limit (τ ′ − τ − 1) → ∞, where A is the amplitude and γ is the dynamical exponent, both to be defined. Supposing that γ < 2 (γ = 2 corresponds to ballistic motion), we expect that both the inner sum (over τ ′ ) and the outer one (over τ ) in eq. (12) will be dominated by the upper summation limit. As a consequence, in the large-t limit and hence, in the large-t limit the expression (12) attains the form In line with the arguments presented in Refs. [27,28], we recall that the dynamical exponent γ defines the characteristic extent of the trajectory Y t ; that being, < Y 2 t >= m 2 t γ , where m 2 is as yet unknown proportionality factor. By symmetry, one expects thus that the DA MSD along the x-axis, i.e., < X 2 t >, obeys exactly the same law, which entails the following closure relation : Inspecting the behaviour of the latter expression in the limit t → ∞, we infer that the contribution of the first term in the rhs of eq. (16) becomes negligible in the limit t → ∞, so that the dominant contribution is provided by the second term. Comparing the power-law on the left-hand-side (lhs) of eq. (16) with the second term on the rhs of this equation, we find that the exponent γ obeys which yields γ = 4/3 -the value which has been previously conjectured and verified numerically in Refs. [27,28]. Therefore, our reasonings correctly reproduce the value of the dynamical exponent γ. However, inferring a numerical value of the prefactor m 2 from eq. (16), (which predicts m 2 ∼ 9A/8), should lead to a somewhat higher m 2 than the actual one, because the rhs in eq. (14) evidently overestimates the value of the double sum in the lhs of this equation. The point is that the algebraic form in eq. (13) is only valid for such realisations of the TP trajectories, for which the sum of the number of jumps and of the number of the pausing events is even. Otherwise, P τ ′ −τ −1 (Y = 0) is exactly equal to zero. As a consequence, eq. (16) overestimates m 2 .
Lastly, we note that a similar type of arguments was invoked to characterise a decay of the number of tree-like clusters with a growing pattern height in a process of ballistic deposition of sticky particles on a line [66]. Both the decay and the ensuing thinning of the forest of such clusters appear to be controlled by a random wandering of the inter-cluster boundaries with the super-diffusive exponent γ = 4/3.
In Fig. 2, panel (a), we present numerical results (open circles) describing the time evolution of the DA MSD of a single TP. The dashed line indicates the superdiffusive power-law behaviour of the form < X 2 t >∼ m 2 t 4/3 , with m 2 = 0.556. This estimate of m 2 is based on the fitting of the full probability distribution, which is discussed below in subsection 3.2. We observe that the super-diffusive behaviour sets in from rather early times and the transient diffusive law, as predicted by the first term in the rhs of eq. (16), is not observed. Next, the inset in the panel (a) illustrates the convergence of the dynamical exponent γ t , defined by to its asymptotic value 4/3. Such a representation of γ t (as compared to the standardly used one, γ t = ln X 2 t / ln(t)) is particularly well-suited for a numerical analysis of the dynamical exponent in an expected power-law dependence on time with an unknown numerical prefactor, since the latter cancels out automatically. In eq. (18) the parameter z is a trial exponent, 0 < z < 1, which rescales time in the second term; in principle, z can be chosen rather arbitrarily; we use z = 0.9. We also observe that γ t converges to its asymptotic value very rapidly, in line with the behaviour of the DA MSD. In panel (b) of Fig. 2 we plot the reduced moments < |X| q t > /t 2q/3 for q = 1, 2, 3 and 4 as functions of time. We observe that the reduced moments saturate as some constant values m q as time progresses, indicating that the moments themselves obey < |X| q t >= m q t 2q/3 (see eq. (20)). Here, the dashed lines indicate our estimates for the values of the numerical prefactors m q (see eq. (21)).

Probability distribution and moments of arbitrary order
In Fig. 3 we depict different facets of the numerically evaluated full probability distribution P t (X, Y ) and of the marginal distribution P t (X), (see eq. (11)). Panel (a) presents the time evolution of P t (X = X ⋆ ) for six fixed values of X ⋆ : X ⋆ = 0, 60, 1000, 1400, 1800 and 3000 (curves from top to bottom, with lighter colours corresponding to smaller values of X ⋆ ). Our numerical results show that, unequivocally, P t (0) obeys a power-law of the form P t (0) ≃ A/t 2/3 , which is fully in line with our above analysis. The decay amplitude is defined with a good accuracy by A ≈ 0.568. Moreover, comparing our numerical results with the form P t (0) ≃ A/t 2/3 , we conclude that the latter provides a very accurate estimate for P t (0) starting from rather short times -the dashed line representing A/t 2/3 and the numerical data (light blue curve) are almost indistinguishable. In turn, P t (X = X ⋆ ) for X ⋆ = 60, 1000, 1400, 1800 and 3000 converges ultimately to P t (0) ≃ A/t 2/3 , which is, of course, not an unexpected behaviour. The panel (b) presents the time evolution of P t (0, 0) -the probability of being at the origin at time moment t. We observe that the power-law form P t (0, 0) ≃ A ′ /t 4/3 (with A ′ ≈ 0.555) describes the numerical data fairly well. Note also that this form implies that a random walk on a random Manhattan lattice is not certain to return to the origin. Further on, in Fig. 3, panels (c) and (d), we plot t 2/3 P t (X) and t 4/3 P t (X, Y ) with Y = 0 as functions of the scaled variable u = X/t 2/3 . The data collapse evidenced by our numerical results for both the central part of the distribution and for its tails, suggests, again rather unequivocally, that the marginal distribution P t (X) of the TP position along the x-axis at (sufficiently large) time t has the following form: where B ≈ 1.249, a ≈ 1.049 and b ≈ 1.730. We observe, as well, that the full distribution P t (X, Y ) (with Y = 0) exhibits essentially the same functional behaviour as a function of u, (see Fig. 3, panel (d)), as the marginal distribution P t (X) and only the values of the parameters are slightly different. We therefore conclude that a) the central part of both distributions is a Gaussian, with the variance which grows super-diffusively with t, and b) the tails of both distributions deviate from a Gaussian and have a form ∼ exp(−|u| 4/3 ), i.e., are "heavier" than a Gaussian. The presence of such tails also manifests itself in the anomalously high asymptotic value ≈ 3.5 attained by the kurtosis of the marginal distribution P t (X) (see the dashed curve in Fig. 7, panel (d)). Recall that the kurtosis of a Gaussian distribution is equal to 3. We note that the large-u tail of P t (X) and P t (X, Y = 0) has a very different form, as compared to the prediction made in Refs. [27,28]. Assuming the validity of the usual relation between the shape exponent δ and the dynamical exponent γ, δ = 1/(1 − γ), it was conjectured that the shape exponent should be δ = 3. Our data shows that this is not the case and, surprisingly enough, the distribution in the second line in eq. (19) has exactly the same shape exponent δ = 4/3 as the one appearing in the MdM model with random layered convection flows (see Refs. [27][28][29]).
Capitalising on the expression in eq. (19), we estimate the behaviour of the moments of P t (X) of arbitrary order q. Multiplying both sides of eq. (19) by |X| q , changing the integration variable for u = X/t 2/3 , and integrating the expression in the first line over u ∈ (−1, 1) and in the second line -over u ∈ (−∞, −1) and (1, ∞), we get with where Γ(a, b) is the incomplete Gamma-function. Note that here we discard the transient region between two asymptotic regimes, supposing that the second regime is valid starting from |u| = 1. This is, of course, not true and hence, m q in eq. (21) overestimates the actual value of the numerical prefactor m q in eq. (20). We however believe that such an estimate is quite plausible and would not incur any significant error. The plot of the numerical results for the first four moments together with the estimates for m q presented in Fig. 2, panel (b), shows that it is indeed the case. We close this subsection with two following remarks: a) the value of m 2 deduced from eq. (16), i.e., m 2 ≈ 9A/8 ≈ 0.639, slightly overestimates the value of m 2 obtained from eq. (21), i.e., m 2 = 0.556. This is completely in line with our argument that the second term in the rhs in eq. (16) provides an upper bound on the actual value of m 2 . b) For the kurtosis κ of the marginal distribution P t (X), i.e., we find from eqs. (20) and (21) that κ = m 4 /(m 2 ) 2 ≈ 3.360, which value favourably agrees with κ ≈ 3.5 deduced from our numerical simulations (see Fig. 7, panel (d)).

Sample-to-sample fluctuations
Up to the present moment we discussed only the averaged behaviour -both over thermal histories and over realisations of random convection flows. However, a legitimate question is how the pertinent parameters themselves vary from a realisation to a realisation of the latter pattern. This question has been first addressed in Refs. [27,28] for the MdM model with random layered flows and significant sample-to-sample fluctuations have been predicted. However, for the dynamics on a random Manhattan lattice this issue has not been analysed and we concentrate on it below, focusing on the mean-squared (averaged over thermal histories only) displacement X 2 t of the TP on a given pattern of arrows as well as on the corresponding probability distribution function Π t (X) of its position X along the x-axis. Recall that P t (X) = Π t (X) .
In Fig. 4, panel (a), we depict the corresponding MSD X 2 t (averaging over thermal histories is performed over 2 × 10 7 realisations of trajectories, for a given realisation of convection flows) for 50 realisations of disorder as functions of time. We do indeed observe some scatter in the values of the prefactor m 2 , which here is a random variable dependent on a particular realisation of disorder. On the other hand, the amplitude of fluctuations does not seem to be very significant and all the curves concentrate essentially around the DA MSD X 2 t = 0.556 t 4/3 (see Fig. 2). Moreover, we realise that the MSD averaged over 50 realisations of disorder only, appears to be fairly close to the DA MSD evaluated using an ample statistical sample; recall that in subsection 3.2 we used 2 × 10 5 realisations of disorder in order to perform averaging over random convection flows. On contrary, the sample-to-sample fluctuations do affect in a significant way the shape of the distribution function Π t (X), both in the central part and especially in the region of anomalous tails, for which the numerical data looks quite nebulous. However, it is quite surprising to realise that being averaged over just 50 realisations of disorder, Π t (X) gets rather close to P t (X) depicted in Fig.  3 (see thin solid and dashed curves in Fig. 4), which again was evaluated using a much bigger statistical sample. Here, an agreement between the blue zigzag curve and the thin black solid line looks nearly perfect within the central, Gaussian part of the distribution, while also for the tails it shows a rather convincing agreement. Therefore, we may conclude that sample-to-sample fluctuations are essentially less important for a random motion on a random Manhattan lattice than in the MdM model with layered random flows [27,28].

Spectral analysis of the TP trajectories
Complementary information about anomalous diffusion of the TP can be inferred from the so-called single-trajectory power spectral density S(T, f ), where T is the observation time and f is the frequency (see, e.g., Ref. [67] for more details). This property is a random, realisation-dependent variable, parametrised by f and T . Note that in the model at hand, it depends on both a given pattern of arrows in a random Manhattan lattice and on a given realisation of a thermal history. For an integer-valued X t , S(T, f ) is defined as and hence, is a periodic function of f T with the prime period 2πT . In a standard text-book analysis, one considers the ensemble-averaged (and also disorder-averaged for our case) value of the random variable in eq. (23), i.e. its first moment : which probes the frequency-dependence of the Fourier-transformed covariance function of X t . Moreover, one also takes formally the limit of an infinitely long observation time, i.e., sets T = ∞. We will demonstrate below that, although one has indeed to consider very large values of T in order to extract a meaningful information about the f -dependence of the power spectral density in eq. (24)), taking the formal limit T = ∞ in our case renders such a standard definition meaningless. In Fig. 5 we present the results of a numerical analysis of the functional form of µ(T, f ), which reveals two rather surprising features. First, it appears that µ(T, f ) is ageing, i.e., its amplitude is dependent on the observation time T , µ(T, f ) ∼ T 1/3 , which is demonstrated in the inset to this figure. As a consequence, setting T = ∞ is meaningless. Second, we observe that f 2 µ(T, f ), (for three fixed values of T ), approaches constant f -independent values for sufficiently large f . This signifies that µ(T, f ) ∼ 1/f 2 , i.e., it exhibits exactly the same f -dependence as the power spectral density of a standard Brownian motion (see, e.g., Ref. [67] and references therein), although the process under study is clearly not a Brownian motion. Therefore, fixing T and focussing only on the f -dependence of µ(T, f ) garnered from numerical simulations, one can be led to an erroneous conclusion that the observed process is a Brownian motion. As an actual fact, this is precisely the T -dependence of µ(T, f ) which helps to realise that this is not the case (see also Ref. [71]).
We note that such a "deceptive" f -dependence has been previously reported for the running maximum of Brownian motion [68], diffusion in a periodic Sinai disorder [69], diffusion with stochastic reset [70] and also for a variety of diffusing diffusivity models [43]. Further on, the law µ(T, f ) ∼ T 1/3 /f 2 was observed for other super-diffusive processes, such as a fractional Brownian motion with the Hurst index H = 2/3 (i.e., γ = 4/3) [71] or a super-diffusive scaled Brownian motion Z t described by the Langevin equationŻ t = t 1/6 ζ t [72], with ζ t being a Gaussian white-noise with zero mean. This latter process also produces a super-diffusive motion with γ = 4/3, suggesting that the law µ(T, f ) ∼ T 1/3 /f 2 might be a generic feature of processes with γ = 4/3. We note parenthetically that this questions the robustness of the textbook approach, based solely on the evaluation of µ(T, f ), which cannot distinguish between these three distinctly different random processes.
The difference between these processes becomes apparent, however, when one considers higher-order moments of S(T, f ), e.g., its variance. In particular, one may focus on the coefficient of variation C v , which is defined by where σ(T, f ) is the standard deviation of a random variable S(T, f ). This characteristic parameter shows a completely different behaviour as a function of f for a super-diffusive fractional Brownian motion and a super-diffusive scaled Brownian motion. For the former C v approaches for sufficiently large T and f a universal (i.e., regardless of the actual value of H > 1/2), time T -independent constant value √ 2 [71], while for the latter -a universal time T -independent constant value √ 5/2 [72], the same which is observed for a standard Brownian motion [67].
To this end, we have studied via numerical simulations the frequency and the observation time dependence of C v for the TP random motion on a random Manhattan lattice. This dependence is presented in Fig. 5, panel (b), in which we plot the coefficient C v of variation of a single-trajectory power spectral density as a function of f for three values of T . We observe that C v tends to a higher than √ 2 value as frequency increases. Moreover, C v is clearly ageing, i.e., its limiting behaviour is dependent on the observation time. In conclusion, we observe a behaviour of C v which is markedly different from the two above mentioned examples of super-diffusion with γ = 4/3.

Tracer particle dynamics on a populated random Manhattan lattice
In this last section we discuss the results of the numerical analysis of the TP dynamics on a crowded Manhattan lattice. In Fig. 6 we present the DA MSD of the TP for Model A and Model B, prefactor m 2 (ρ) in the super-diffusive law X 2 t = m 2 (ρ) t 4/3 as a function of the density ρ of the LG particles, and also the kurtosis of the distribution P t (X) for Model A and Model B. First of all, we realise that for both models the DA MSD obeys the same super-diffusive law X 2 t ∼ t 4/3 , for any density of the LG particles. Prefactor m 2 (ρ) depends, of course, on the density of the LG particles and their dynamics; indeed, m 2 (ρ) shows apparently different dependences on ρ for Model A and Model B. For Model A, in which the TP follows random convection flows while the LG particle perform constrained random walks, this dependence is most strong and is rather close to a parabolic law m 2 (ρ) = 0.556 (1 − ρ) 2 (thin solid curve in panel (c)). This parabolic law very accurately describes the actual dependence of m 2 (ρ) on ρ for ρ < 1/2. For higher densities, however, some deviations are clearly seen, although such a discrepancy can be also attributed to the lack of a large enough statistical sample. For Model B, in which all particles are identical and all perform a super-diffusive motion, the ρ-dependence of the prefactor m 2 (ρ) is given by m 2 (ρ) = 0.556 (1 − ρ) 4/3 (thin dashed curve in panel (c)). This law agrees fairly well with the numerical data for any value of the LG particles density. In turn, as we have mentioned above, such a dependence implies that the system is perfectly stirred and the time variable t gets merely rescaled by the fraction of successful jump events of the TP, i.e., by (1 − ρ), as one can expect from simple mean-field-type arguments. Arguably, this expression for m 2 (ρ) is exact.
Next, in Fig. 6, panel (d), we depict the time evolution of the kurtosis of the distribution P t (X) in case of a single TP, as well as for Model A and Model B. Thick dashed line corresponds to the kurtosis of P t (X) in case of a single TP. We observe that κ saturates as t evolves at a constant value which is close to 3.5, i.e., it exceeds the value 3 specific to a Gaussian distribution and thus implies that P t (X) is not Gaussian. This happens, of course, due to the presence of anomalous slower-than-Gaussian tails, which we discussed in subsection 3.2. Next, noisy green curves depict our results for κ in Model A, with the LG particles densities ρ = 0.1, 0.2, 0.3 and ρ = 0.4, plotted versus a rescaled time variable t ρ = (1 − ρ) 3/2 t. Although there is a significant scatter of these curves at short times, we notice that for sufficiently large t ρ all these curve collapse on the dashed line representing a single TP case. For Model B we analyse the TP dynamics for a broader range of the LG particles densities; we considered nine values of ρ, ρ = 0.1, 0.2, . . . , 0.9. Here, the curves defining the evolution of κ for different values of ρ, merge altogether even for short times when plotted versus a rescaled time t ρ = (1 − ρ)t, and eventually approach the value of the kurtosis for a single TP case. Such a behaviour signifies that for sufficiently large times P t (X) possesses some universal scaling properties for both Model A and Model B.
Lastly, in Fig. 7 we present the numerical data for the distribution function P t (X) for Model A and Model B. We observe that for both models the central Gaussian part of the distribution is described with a very good accuracy by a Gaussian function: where the choice of t ρ depends on the model under study. For a single TP case, t ρ = t. The lack of a big statistical sample does not permit to make completely conclusive statements about the tails of the distribution. We notice, however, that such tails are definitely present and a departure of the distribution from a purely Gaussian form is apparent in Fig. 7. We also observe that upon an increase of t the curves get closer for both Model A and Model B. We thus find it absolutely plausible that such a form of distribution is also valid for the dynamics of a TP on a populated random Manhattan lattice.

Conclusions
To recapitulate, we studied the tracer particle (TP) dynamics in presence of two interspersed and competing types of disorderquenched random convection flows on a random Manhattan lattice, which prompt the TP to move super-diffusively, and a crowded dynamical environment formed by a lattice gas (LG) of hard-core particles, which hinder the TP motion. The random Manhattan lattice is a square lattice decorated with arrows in such a way that directionality of each arrow is fixed along each raw (a street) or a column (an avenue) along their entire length, but whose orientation randomly fluctuates from a street to a street and from an avenue to an avenue.
The hard-core LG particles perform a random motion, constrained by the singleoccupancy condition; that being, each lattice site can be occupied by at most a single particle -a LG particle or a TP, or be vacant. We have considered two possible scenarios of the LG particles random motion. In Model A, we supposed that the LG particles are insensitive to the random convection flows and perform symmetric random walks -a simple exclusion process -among the sites of a two-dimensional square lattice. In this case, the TP moves subject to random convection flows and interacts with a fluid-like quiescent environment, which imposes some frictional force on it. In Model B, we supposed that all the particles -the TP and the LG particlesare identical and follow a local directionality of bonds in a random Manhattan lattice. In this case, the system under study is a kind of a "turbulent" fluid in which all the particles perform a super-diffusive motion.
We focused on such characteristics of the TP dynamics as its disorder-averaged mean-squared displacement (DA MSD), (and generally, the moments of arbitrary order for a single TP), the distribution of its position at time moment t averaged over disorder, and the time evolution of the kurtosis of this distribution. We have shown that for both Model A and Model B the DA MSD obeys a super-diffusive law X 2 t ∼ m 2 (ρ) t 4/3 , where the functional dependence of the prefactor m 2 (ρ) on the mean density ρ of the LG particles depends on the model under study. For the case of a single TP (i.e. in absence of the LG particles, ρ = 0), we provided some analytical arguments explaining such a super-diffusive behaviour.
We showed that the distribution of the TP position has a Gaussian central part and exhibits slower-than-Gaussian tails of the form exp(−(|X|/t 2/3 ) 4/3 ) for sufficiently large X and t. Such a form was evidenced in case of a single TP through an analysis of a very big statistical sample, and also shown to hold, although not in a completely conclusive way, for the dynamics on a populated random Manhattan lattice. As a consequence of presence of anomalous tails, the kurtosis of the distribution in all the situations under study, was shown to attain a bigger value (3.5) than the value (3) specific to a Gaussian distribution.
Finally, we addressed the question of sample-to-sample fluctuations in the system under study and performed an analysis of spectral properties of the TP trajectories, which revealed some interesting features.