Diffusion in time-dependent random environments: mass fluctuations and scaling properties

A mass-ejection model in a time-dependent random environment with both temporal and spatial correlations is introduced. When the environment has a finite correlation length, individual particle trajectories are found to diffuse at large times with a displacement distribution that approaches a Gaussian. The collective dynamics of diffusing particles reaches a statistically stationary state, which is characterized in terms of a fluctuating mass density field. The probability distribution of density is studied numerically for both smooth and non-smooth scale-invariant random environments. Competition between trapping in the regions where the ejection rate of the environment vanishes and mixing due to its temporal dependence leads to large fluctuations of mass. These mechanisms are found to result in the presence of intermediate power-law tails in the probability distribution of the mass density. For spatially differentiable environments, the exponent of the right tail is shown to be universal and equal to −3/2. However, at small values, it is found to depend on the environment. Finally, spatial scaling properties of the mass distribution are investigated. The distribution of the coarse-grained density is shown to possess some rescaling properties that depend on the scale, the amplitude of the ejection rate and the Hölder exponent of the environment.


Introduction
Many different situations found in nature occur in an environment where spatial fluctuations occur at length scales so small that they contribute to macroscopic processes only through an averaged effect. The environment is then usually modeled by introducing some random disorder from which effective properties are derived. Examples include the Kraichnan velocity ensemble for turbulent transport [1] and the Sherrington-Kirkpatrick spin-glass model for magnetization [2]. There are many applications where the environment fluctuates on much slower timescales than the processes of interest (such as diffusion, transport or wave propagation). This leads us to consider quenched disorder, as in the case of directed polymers [3,4], of wave propagation in random media [5], of heat transport in open-cell foams [6] and of viscous permeability in porous materials [7]. In most instances, one of the main goals is to derive an effective diffusivity, conductivity or permeability of the medium (see, e.g., [8]).
As to diffusive processes in quenched random media, mathematicians and physicists have largely studied their discrete version known as random walks in random environment (RWRE) [9]. One usually considers a lattice with fixed random transition probabilities between sites and studies the behavior of a random walk on it. One of the important questions of interest has been to determine under what assumptions such RWRE are transient or recurrent (i.e. whether they escape to infinity or indefinitely come back to their starting position). For a time-independent random environment with no space correlations, it has been shown that randomizing the environment slows down the diffusive properties of the random walk [10]. Furthermore, under some precise assumptions, Sinai [11] proved that one-dimensional (1D) symmetric nearest-neighbor random walks in an uncorrelated environment are sub-diffusive and scale as X t ∼ log 2 t when t → ∞. Derrida and Pomeau [12] demonstrated that when the assumption of symmetry is relaxed, the random walk escapes as X t ∼ t α , with 0 < α < 1. Bricmont and Kupiainen [13] showed that the upper critical dimension of random walks in a time-independent random environment is 2, so that for d > 2 RWRE are diffusive. More recently, a great deal of effort was devoted to proving central-limit theorems or large-deviation principles for quenched disorder [14][15][16]. Such works generally assume that the environment has no (or very short) spatial correlations and that it is time independent. Very little is known in 3 the time-dependent environment. One then deals with annealed statistics where the averages are performed over the fluctuations of the random medium. In such settings and under some rather general assumptions, it was shown that random walks are always diffusive and that central-limit theorems generally hold [17,18].
In this work, we are interested in those situations where the timescales and length scales of diffusion are comparable to or longer than those at which the environment fluctuates, and thus we are interested in time-dependent random media with both temporal and spatial correlations. A clear instance where this occurs is turbulent transport whose modeling is a key issue in industrial and environmental science. The classical models generally used in engineering and meteorology are based on an 'eddy-diffusivity' approach (see, e.g., [19]), which in general requires a clear scale separation between the turbulence and large-scale variations of the mean velocity. Homogenization techniques can then be used to show that the averaged concentration field follows an effective advection-diffusion equation with an advecting velocity and a diffusion coefficient tensor that depend on the slow variables [20]. The large-scale mixing originating from small-scale fluctuations acts only through the diffusive term, which, by the maximum principle, cannot be responsible for the creation of large concentrations. The concentrations observed in compressible flows can thus come only from the effective advection term. As we will see later in this paper, the situation can be rather different when one is interested in other types of compressible transport.
The paper is organized as follows. In section 2, we introduce a time and space continuous model of diffusion in a random environment that can be interpreted either as a limiting process of Langevin dynamics or as the continuous limit of a discrete ejection model. The model corresponds to an Itō differential equation with multiplicative noise and no drift. Section 3 is devoted to the study of individual trajectories of the Itō diffusion in the case of smooth and non-smooth environments. When the environment has a finite correlation length, particles are found to have standard diffusion properties at large times with a displacement distribution that approaches a Gaussian. In section 4, we define the density of mass and study the fluctuations of the density in smooth and non-smooth random environments. The locations where the environment ejection rate vanishes are shown to play a crucial role in the large fluctuations of density. Spatial scaling properties of the mass distribution are then investigated in section 5. We show that the distribution of the coarse-grained density possesses some rescaling properties that depend on the scale, the amplitude of the ejection rate and the Hölder exponent of the environment. Finally, section 6 presents some concluding remarks.

Description of the model
We consider the dynamics given by the following Itō stochastic differential equation: where W t is the D-dimensional Wiener process. The diffusion coefficient σ is a random, spaceand time-dependent field, whose statistics are independent of W t and are described below. Before specifying our settings, we motivate this model by considering two instances where (1) arises as a limiting process.

4
Let us first consider species whose dynamics is dissipative in the full position-velocity phase space. We assume that the trajectory of a particle obeys the Newton equatioṅ The forces acting on the particle are viscous drag and an external force f that depends on both space and two different timescales. We have dropped here the vectorial notation but all the following considerations can be easily generalized to any dimension. Such an equation describes for instance the dynamics of a heavy inertial particle in a velocity field that varies over two timescales. Also, one could consider a large particle embedded in a time-and space-dependent thermal bath. The fast timescale ε can be interpreted as the typical time of momentum exchange between the particle and its environment. As in Einstein's original work on Brownian motion, the timescale ε is much smaller than the viscous damping time 1/µ. In the limit εµ → 0, the force f can be approximated by a Gaussian noise with correlation where the average · ε is with respect to the fast time variable. The spatio-temporal variations of σ 2 can be interpreted as a non-homogeneous temperature field in the thermal bath. We next consider the limit when the response time 1/µ is much shorter than the slow timescale. This introduces a new fast timescale O(1/µ) over which the particle velocity fluctuates. We see that where t ∧ t = min(t, t ). The last equality is obtained by integrating by parts and neglecting exponentially small terms. Then, taking the limit µ → ∞ in (4) we obtain V t V t = 2 σ 2 (X t , t)δ(t − t ). The particle position thus satisfies the stochastic equation (1). The model we consider is thus relevant to the study of the Langevin dynamics in a fluctuating environment. The second instance where (1) arises is in the continuous limit of a mass-ejection model. Let us consider the D-dimensional periodical lattice of fundamental size x and period L = N x. This lattice defines a tiling on which we consider the following discrete-time dynamics. The cell indexed by i contains at time n t a mass ρ i (n). Then, between times n t and (n + 1) t, a fraction 0 < γ i (n) < 1/(2D) of this mass is ejected from the cell i to one of its 2D neighbors. For a given set of ejection rates {γ i (n)} i∈{1,...,N } D the variables {ρ i (n)} i∈{1,...,N } D define a Markov chain whose master equation is given by where N i are the neighbor cells of i. It is clear that, because of periodic boundary conditions, the total mass M = i ρ i is conserved. Equation (5) describes a mass-ejection process in a timedependent non-homogeneous environment determined by the ejection rates γ i (n). A similar model has been studied in [21] in the case when the γ i (n) take two values, γ and 0, and are uncorrelated in both space and time. This model was motivated by the study of inertial particles clustering in turbulent flows and was expected to mimic heavy particle ejection from coherent vortical structures by centrifugal forces.
Turning back to the ejection model (5), we now consider the continuous limit x → 0 and t → 0. To take this limit, we suppose that the ratio x 2 / t is kept fixed. We denote by σ 2 (x, t) the continuous limit of γ x/ x (t/ t), where the coefficient is the typical amplitude of the ejection rate and σ is an order-one time and space continuous function. In this limit, equation (5) becomes This is the Fokker-Planck (forward Kolmogorov) equation associated with the Itō stochastic differential equation (1). For t > s, we define the (forward) transition probability density where · designates averages with respect to the realizations of W t . Then the solutions to (6) trivially satisfy for all t > s The model described by the diffusion equation (6) can thus be reinterpreted in terms of an RWRE. Indeed the discrete version of the process given by (1) corresponds to a random walk on the D-dimensional lattice, with a probability to hop from site i to one of its neighbors equal to 2D γ i (n).
In contrast with a standard diffusion equation coming from Fick's first law, the differential operator ∇ 2 [σ 2 (x, t) ·] in the right-hand side of equation (6) is not positive-definite. There is no maximum principle and the solution is expected to behave very differently. The mass is going to accumulate at the zeros of σ . Suppose indeed in one dimension that σ (x, t) C x in the vicinity of x = 0. At leading order, the flux reads J = − (∂/∂ x)(σ 2 ρ) −2 C 2 xρ, which is positive for x < 0 and negative for x > 0, leading to a permanent mass flux toward x = 0. As we will see later, the zeros, their densities and their lifetimes play a crucial role in the statistical properties of the density field.
We now specify how the statistical properties of the random environment are prescribed in our model. We consider that L = 2π and that σ (x, t) is 2π -periodic in space. The environment is then entirely determined in terms of its Fourier series that we write as where the a k are positive real amplitudes and the τ k are scale-dependent characteristic times. The Fourier modes χ k satisfy χ −k (t) = χ k (t) and χ 0 ≡ 0. They are independent Gaussian processes with unit variance and unit correlation time. We choose them as complex Ornstein-Uhlenbeck processes that solve the stochastic differential equation [22] dχ k (s) = −χ k ds + √ 2 dB k (s), where the B k are independent 1D complex Wiener processes. We now prescribe scale-invariance properties for the random environment. We want |σ (x + r, t) − σ (x, t)| ∼ |r| h , where h > 0 controls the regularity of the ejection rate. When h < 1, it is the spatial Hölder exponent of the ejection rate. This amounts to assuming that the Fourier mode amplitudes behave as a k = 1 2 k −1/2−h . Temporal scale invariance is set by assuming that τ k = k −β with β > 0. The exponent β relates time dependence to spatial scale invariance; the case β = 0 corresponds to σ having a unique correlation time and β → ∞ to white noise. As we shall see below, when h < 1, dimensional analysis motivates the choice β = 2 − 2h.
The numerical results presented in this work are obtained in two different ways. The density evolution is obtained by solving equation (6) with a second-order finite-difference discretization of the Laplacian and a semi-implicit Euler temporal scheme. This numerical scheme ensures the conservation of the total mass with high accuracy. In this work, numerical resolutions vary from N = 512 to N = 8192 collocation points and the time step is chosen small enough to well resolve the fastest timescale of the problem. For particle trajectories, the stochastic equations (1) and (10) are solved by a standard stochastic Euler scheme and ensemble averages are obtained by Monte Carlo methods. In both cases, we expect the error made in the solution to act as a numerical diffusion with a constant proportional to the time step. Figure 1(a) show the temporal evolution of the density and of the ejection rate σ 2 (x, t), respectively, for a typical realization in dimension D = 1. Figure 1(b) shows a snapshot of the random environment σ and of the position of a set of particles obeying equation (1) in D = 2. It is apparent in figure 1 that the high-density zones are located in the vicinity of the zeros of σ (x, t). From the right of figure 1(a) we observe that the zeros of σ follow random paths. Their role in the statistical properties of the density of mass and in the diffusion properties of particles will be discussed further in this work.
The distribution and time evolution of these zeros strongly depend on the parameter h characterizing the environment. Let us first consider the case of large h. In figure 1(a), it is apparent that the zeros typically appear in pairs and then diffuse and separate until they merge together or with another zero. Without loss of generality, we can assume that the ejection rate 7 has only two modes and takes the simple form with σ r + i σ i = 2χ 1 . Using the Itō formula it is possible to show that σ can be equivalently rewritten as σ ( where the stochastic processes A t and t satisfy where B A t and B t are uncorrelated Wiener processes. It is clear that the amplitude A t fluctuates around A = 1 and that its correlation time is of the order of 1. It follows from equation (12) that, for typical values of A t , the phase t , and thus the positions of the zeros of σ , diffuse on a timescale of the order of 1.
When h < 1, the random environment presents scaling properties and very different behavior is expected. We have by construction that the second-order structure function of σ behaves like where the overline stands for the ensemble average with respect to fluctuations of the environment (i.e. with respect to the Ornstein-Uhlenbeck processes χ k ). In this work, we make the choice of relating space and time correlations by a dimensional argument. The characteristic time τ k introduced in equation (9) behaves as a power law, so that the correlation time of σ at scale is τ C ( ) = ( /L) β . According to equation (13), the diffusion timescale associated with the spatial scale behaves as τ β+2h−2 defines a dimensionless space-dependent parameter that is usually called the Kubo number. When Ku( ) 1, the environment looks as if frozen. When Ku( ) 1, it fluctuates in an almost time-uncorrelated manner. When β = 2 − 2h, the Kubo number depends on the scale and this breaks any possible scale invariance of the mass distribution. When β = 2 − 2h, we have Ku( ) = Ku = /L 2 , so that scale invariance is possible. Here, we will focus on the latter case and describe the mass concentration properties in terms of the dimensionless parameter Ku ∝ . Note that the choice of having a scaleindependent Kubo number is common in the framework of turbulence and is expected to be relevant to the problem of diffusion of inertial particles in turbulent flows.
Finally, let us make a couple of remarks on the dependence upon another parameter of the environment, namely the Hölder exponent h. Using Parseval's theorem, it is possible to show that σ 2 is a Gaussian variable with the variance of increments given by equation (13). Therefore σ (·, t) is a fractional Brownian motion of exponent h [23]. For h = 1/2 it corresponds to standard Brownian motion. For h > 1/2 the increments of σ are not independent and have positive covariance; the zeros of σ are then finite and isolated. For h < 1/2, the covariance of increments is negative, and therefore we expect the number of zeros to become infinite and to accumulate. The different behavior of the random environment will affect the general properties of diffusion, as we will see in the following sections.
As the Ornstein-Uhlenbeck processes χ k are stationary, the fluctuations of the ejection rate σ 2 are a stationary and homogeneous random field. Hence we expect that, at sufficiently large times, the density distribution reaches a statistically stationary state and numerical simulations 8 indicate that this is indeed the case. The work reported in this paper is mainly concerned with the statistical properties of the density field in this large-time asymptotics. Most of the results presented in this work are in dimension D = 1.

Diffusive properties
We now turn to studying the diffusive properties of the solutions X (t) to (1). One has so that where the angular brackets denote average over trajectories. At times much larger than the correlation time of σ 2 along particle trajectories, the integral becomes a sum of independent random variables, which obeys the law of large numbers. We thus have where the overline denotes the average with respect to the environment. The effective diffusion constant D( , h) involves an average of the environment along particle trajectories which is equivalent to an average weighted by the mass density ρ. As the trajectories are expected to spend a long time in the vicinity of the zeros of σ , we expect this Lagrangian average to be less than the full Eulerian spatial average. The displacement fluctuations of order less than √ t are expected to be given by a central-limit theorem. However, larger fluctuations should obey a large deviations principle with a rate function that might not be purely quadratic as it is nontrivially related to the Lagrangian properties of the environment. To study these fluctuations, we perform Monte-Carlo simulations of equation (1).
We first consider the case when the ejection rate σ 2 is a smooth function of space (i.e. h > 1). With this configuration there is a finite (and small) number of isolated zeros separated by a distance order L/2 = π . A number of simulations of equation (1) have been performed for different values of using 32 modes and h = 2. Figure 2 shows the meansquare displacement of particles averaged over 1000 realizations of the diffusion and 20 000 realizations of the environment. As expected from equation (16), linear growth is observed at large times. However, the diffusion constant D( ) non-trivially depends on , as seen in the inset of figure 2. The non-monotonic behavior of D( )/ can be explained by two competing phenomena. In the limit → 0 (or equivalently Ku → 0) the environment changes fast enough compared to the characteristic time of diffusion given by 1/ . In this limit, we can consider that the particle trajectories sample the environment homogeneously, so that the Lagrangian average can be replaced by the spatial Eulerian average. This leads to  (16) on . For h = 2, we obtain D( , 2)/ 1.03, which is in agreement with the value observed in the inset of figure 2. As increases, trajectories spend a longer time at the zeros of σ decreasing the value of D( )/ . In the limit of → ∞ (i.e. Ku → ∞), the environment can be considered completely frozen and we expect the mass density solving equation (6) to be approximately given by ρ(x, t) ∼ 1/σ 2 (x, t). From equation (16), we obtain that D( , h) = ( /L) σ 2 ρ d(x) ≈ ( /L) dx = . Hence D( )/ = 1 for → ∞ as seen in the figure. In this limit the mass is completely concentrated at the zeros of the ejection rate, and the diffusion is carried out by the motions of the latter that follow a dynamics similar to equation (12).
The normalized probability distribution functions (PDFs) of the displacements are shown in figure 3 for different times and different values of . At early times and for low values of , the PDFs present exponential tails that are also observed at intermediate times for larger values of . In all cases, the PDFs approach a Gaussian distribution at very large times. Note that for sufficiently large values of , some oscillations can be observed on the PDF tails at intermediate times. As we will now see, this is due to trapping by the zeros of σ (x, t). To emphasize this point, we compare the PDF of the displacement for different values of chosen at the time t * such that (X (t * ) − X (0)) 2 = (L/2) 2 = π 2 , i.e. the time when the typical distance traveled by the particles reaches a length of the order of the separation between two zeros (of the order of L/2 for smooth environments). The PDF of (X (t * ) − X (0))/(L/2) is displayed in figure 4(a) for various values of . It is apparent that the bumps appear at multiples of L/2 for large s. With a high probability, a particle is initially located close to a zero of σ . When it travels away, there is a very good chance that it stops at the neighboring zero, leading to a quasi-multimodal distribution.
To better quantify the large-time behavior of the displacement, we compute its higher-order moments and compare them to those corresponding to a Gaussian distribution. We thus define for even p where [z] denotes the Gamma function. C p = 0 for all p 2 if X (t) − X (0) is a Gaussian variable. C 2 = 0 by construction and C 4 is the kurtosis. C p can be interpreted as a pth-order deviation from a Gaussian distribution. Data suggest that all the C p display 1/t behavior at large times. The dependence on the amplitude and on the order p is well represented by the functional form The temporal evolution of C p /A( ) for different values of and p, and the function A( ) are displayed in figures 4(b) and (c), respectively. As can be seen from figure 4(b), the data shows the dependence of the diffusion constant normalized by D 0 = σ 2 L 2 = ζ (1 + 2h). It is seen that the importance of zeros on the diffusion constant increases as h decreases.

Density fluctuations
So far, we have studied the diffusion properties of particles described by the stochastic equation (1). We now turn to studying the statistical properties of the mass density ρ that evolves according to the diffusion-like equation (6). As we have previously observed, the spatial fluctuations of ρ are correlated with the distribution of zeros of the random environment σ (x, t). More precisely, we observe from equation (6) that the flux of mass at x is directly given by The flux hence vanishes at the zeros of σ (x, t), so that they would act like sinks concentrating all of the mass as if they were not moving. However, due to their diffusion, the zeros are not able to indefinitely concentrate mass and the density saturates in their neighborhood. In parallel, as the total mass is conserved, this concentration process leads to the creation of voids in the regions where σ is of the order of 1.
In the following, we consider the case D = 1. Under some assumptions of homogeneity and ergodicity of the dynamics, the PDF of the mass density can be written as This will be the main quantity studied in the following sections. As we will see, it strongly depends on the parameters and h. The definition (22) in terms of a space-time average will be particularly useful in first attacking the stationary case and then developing phenomenological arguments on its behavior in the non-stationary case.

Smooth random environments
Let us first consider the limit Ku → ∞, i.e. when the ejection is infinitely larger than the inverse of the diffusion time of the zeros of σ . Note that, in the limit of a time-independent environment, no stationary state can be achieved. Mass will then concentrate around the zeros, and at t = ∞ the density will become atomic, supported at these points. For sufficiently large but finite times, as almost all the mass is concentrated around the zeros of σ , the ejection rate can be Taylor-expanded in the vicinity of these points. Without loss of generality, let us assume that a zero appears at the origin x = 0 at t = 0 and does not move. In this limit, equation (6) reduces to where C = dσ/dx| x=0 . Rescaling time allows us to set C 2 = 1. Without loss of generality, we also consider an infinite domain and an initial condition with compact support: ρ(x, 0) = 1 for |x| 1 and 0 elsewhere. Equation (23) can be analytically integrated by using the method of characteristics. For x > 0 the change of variable u(y, t) = e 2t ρ(x = e y−3t , t) leads to the homogeneous heat equation that can be easily solved by using the corresponding Green function. The solution for x < 0 is obtained by symmetry. One then obtains where erfc(z) = 2 √ π ∞ z e −s 2 ds is the complementary error function. The time-averaged distribution P T (ρ) defined in equation (22) can be rewritten as where xρ(t) is such that ρ(xρ(t), t) =ρ. It exists only for t > lnρ/2 . The contour lines of ρ(x, t) are displayed in figure 6(a). The high-density zones are concentrated in a narrow region of space. This will justify the saddle-point approximation used below. With the change of variables t → λ = 2t/ lnρ, introducing µ(λ;ρ) = − ln xρ(t)/ lnρ, and after some algebra, equation (25) becomes with F(λ;ρ) = λ + µ(λ;ρ) − 1 2λ where in the last equation the asymptotic of erfc −1 (z) ≈ log (1/z) for z ≈ 0 has been used (recall that λ > 1). Using a saddle-point approximation for lnρ 1 of the integral in equation (26), we obtain at leading order (dropping the bar over ρ) This scaling is clearly observed in figure 6(b), where equation (6) is integrated using the stationary environment σ (x, t) = sin x.
Note that a clear power-law scaling is also seen in figure 6(b) at small values of the mass density. To explain this scaling, let us consider a cell of length x located at x 0 far away from a zero of σ (x) and where the density reaches its minimum. At leading order the mass ejected from this cell between t and t + dt is proportional to J ρ(x, t) ∂ x σ 2 (x)| x=x 0 as ∂ x ρ ≈ 0 at this point. Considering that ∂ x σ 2 (x) is constant near x 0 yields an exponential decay of the mass in the cell. Introducing in equation (22) a density of mass with an exponential decay ρ(x, t) ∼ ρ 0 e −σ 2 (with σ 2 0 an effective rate) directly leads to which is in agreement with figure 6(b). As we will soon see, the situation is rather different in the non-stationary case where the motion of the zeros of the ejection rate limits the process of mass concentration. Now, the coefficients σ r (t) and σ i (t) of equation (11) are the Ornstein-Uhlenbeck processes defined in equations (9) and (10). The PDF of ρ for different values of is presented in figure 7. The time needed for the density to accumulate near the zero of the ejection rate σ 2 is of the order of −1 .
Hence, for small values of (i.e. Ku 1), the fast temporal variations of σ (x, t) do not leave enough time for the system to accumulate mass. It is apparent in this case that most of the mass is distributed around the mean mass ρ = 1 (see the = 0.1 curve represented by orange stars). However, for large but finite values of Ku, accumulation is fast enough and a resulting ρ −3/2 scaling is observed. This scaling can be derived by considering that for large Ku the system rapidly relaxes to a quasi-stationary solution in a time of the order of 1/ and stays there for a time of the order of the correlation time of σ . Far enough from a zero, the environment behaves like σ (x) C x and its time evolution can be neglected. The density converges there to a quasistationary solution, which is such that ∂ 2 x (σ 2 ρ) ≈ 0, so that ρ ≈ σ −2 ∝ 1/x 2 . Closer to a zero, the time variations of its location become faster than the concentration of mass and the density saturates to a finite value ρ max (see the left-hand side of figure 8). The transition between these two regimes occurs at a distance x from a zero, which by continuity satisfies x ∼ ρ −1/2 max . The length x is of the order of the distance traveled by a zero during a timescale equal to that of mass concentration. Hence, if we assume that a zero diffuses, one has x ∼ −1/2 , so that the typical value of ρ max is of the order of . The distribution of density has thus a crossover at ρ ∼ . When ρ , that is, for values much below the plateau at ρ max , the behavior is dominated by the divergence of the quasi-stationary density profile. Introducing ρ(x) ∼ 1/x 2 in equation (22) straightforwardly leads to the algebraic behavior 2 max x x 1/ Figure 8. Left: sketch of the density field in the vicinity of a zero of the ejection rate σ 2 located at x = 0; the density grows as x −2 and saturates to ρ max at a distance x ∼ ρ −1/2 max from the zero. Right: sketch of the diffusive time evolution of a zero of σ 2 , which has remained at a distance less than x for a time T before the reference time t = 0.
Note that this value of the exponent in the large-density intermediate tail is very robust. It has also been observed for larger numbers of modes and with an environment of the form σ (x, t) = √ cos (x − ct) (data not shown). The behavior at ρ of the PDF of density is related to the large fluctuations of ρ max and thus to the events when a zero does not move much during a time greater than −1 . More precisely, if at a fixed time we observe a large value of ρ max , this means that the zero has not moved by a distance larger than ρ −1/2 max during a time interval of length T larger than −1 (see the right-hand side of figure 8). When zeros diffuse, this probability relates to the first exit time distribution of a diffusive process. The probability density of the first time T at which the Wiener process exits the interval |x| < x is ∼ exp(−C T / x 2 ) (see, e.g., [24]). We thus obtain For ρ the leading-order behavior is given by T = −1 . This leads to the following exponential cutoff for the PDF of mass density: This exponential behavior is confirmed numerically as can been seen from the lin-log plot in the inset of figure 7.

Non-smooth random environment
We now turn to the case of a time-dependent non-smooth ejection rate σ . Figure 9 presents the PDF of ρ for a number of runs with different values of and h. Note that for large power-law behavior is observed at small masses-see the inset of figure 9(b). These tails can be understood using a phenomenological argument similar to the one used in [21]. Consider an extreme event leading to a very low density in a cell where only ejection of mass has taken place for some time.
For such a configuration we expect an exponential decay of the density and thus a time of the order of T ∝ − −1 ln ρ to reach a small density ρ (where σ 0 is an effective rate). The probability of having this configuration depends only on the properties of the environment. Let p < 1 be the probability of having such a configuration. As the environment de-correlates in a time of the order of 1, the number n of independent realizations of the environment that is needed for the full process to take place during a time T is proportional to T , and thus n = −C −1 ln ρ, where C is a constant that depends on h. Therefore, the probability of this complete event is We thus obtain that the probability of such an extreme event is ρ β , with β = −C ln p/ > 0. Note that this is actually a lower bound to the small mass tail, which is valid as far as p = 0 and C = 0. Obtaining an analytic expression for the exponent β is a theoretically challenging problem that is beyond the scope of this work. The exponent β can be obtained by fitting the left tail of the mass density PDF inside the range where the power law is observed. The fits are presented as dashed lines in the inset of figure 9(b), and the dependence of β on the Hölder exponent h is displayed in figure 10(a) for = 10.
Note that the width of the PDFs in figure 9 clearly depends on the value of the Hölder exponent h as is apparent in figure 10(b), which represents the variance of the density δρ 2 = (ρ − ρ ) 2 . Note that neither the exponent β nor the variance δρ have monotonic behavior as a function of h. This can be qualitatively understood in the following way. When the value of h is decreased, the number of zeros of σ (x, t) increases. However, their spatial distribution also depends on the Hölder exponent h. On the one hand, for h > 1/2 the space increments of σ are positively correlated. This implies that there is no accumulation of zeros. For this configuration, the mass is transported from the non-vanishing zones of σ (x, t) to the nearest zeros. Therefore as h decreases, more and more zeros are present, creating large masses and void zones and thus increasing fluctuations. On the other hand, for h < 1/2 the covariance of increments is negative and there are some finite-size regions with a large number of zeros. σ (x, t) vanishes many times in such zones, so that the diffusion is very weak and mass is trapped. This makes the transfer inefficient and reduces the probability of having extremely high or low mass concentrations. In simpler terms, a large number of zeros increases fluctuations as long as they are not too dense. We expect then a change of behavior near h = 1/2; this is in agreement with figures 10(a) and (b), especially for the largest value of where the properties of the system are expected to depend more strongly on the ejection rate distribution of zeros.

Scale invariance of the mass density field
In this section we study the scaling properties of the density field. For this, we introduce the coarse-grained mass density Note that by homogeneity we have ρ = ρ . The PDF of ρ is computed using definition (22) and is displayed in figure 11 for different values of h and = 1. The narrowest curves on each figure correspond to the largest (non-trivial) scale L/2 and the curves monotonically increase in width with decreasing scale, down to the smallest scale of the runs. For spatially smooth environments, the coarse-grained density of mass becomes invariant when decreasing the scale, as is apparent in figure 11(a) for h = 2. This indicates that the density field is spatially smooth and that the zeros are isolated. The collapse is faster for small densities. This is due to the accumulation of mass only in some very small clusters and the creation of large void zones. These large voids dominate the coarse-grained density statistics up to their typical size, that is, for < L/8, as seen from figure 11(a). In contrast, the typical size clust of mass clusters is rather small. Averaging over scales larger than that, this length scale reduces the largest mass fluctuations. Indeed, the coarse-grained density cannot exceed clust / and this cutoff decreases with , as can be observed in figure 11(a).
This argument does not apply for non-smooth environments where the distribution of zeros is highly non-trivial. The scale invariance is then broken, as can be observed in figures 11(a) and (c). This is also apparent in figure 12(a) figure). However, one observes that for h < 1 and L the variance presents a power-law scaling δρ 2 ∼ −ζ that suggests self-similar behavior of the density. The exponent ζ (h; ) is displayed in the inset of figure 11(c).
The scaling invariance of the density field can be understood in the neighborhood of the zeros of σ . Let us assume that the ejection rate vanishes at x = x 0 . Then, in the neighborhood of x 0 , we have σ (x, t) 2 = |σ (x 0 + x, t)| 2 ∼ x 2h . Supposing that the density behaves as a power-law in the vicinity of x 0 and rescaling space asx = λx, one can easily see that the rescaled density Hence, we expect that if the coarse-grained density ρ presents a self-similar property, then the PDF of ρ obtained with a given value of will coincide at large values with that of ρ λ corresponding to an ejection rate of amplitude λ 2−2h . This scale invariance is apparent in figure 12(b) where different PDFs of ρ , with different values of and such that 2h−2 = const, are confronted for h = 1/2 and collapse at large densities.

Conclusions
To summarize, we have introduced a diffusion model based on a simple discrete mass ejection process where the ejection rates are random variables with temporal and spatial correlations. The model can be interpreted in terms of random walks in a time-dependent random environment. We have considered space-periodic environments where the temporal dependence of each Fourier mode is given by independent Ornstein-Uhlenbeck processes and both the amplitudes of the modes and their correlation times present some prescribed spatial scaling properties. This allowed us to consider smooth and non-smooth environments with fast and slow temporal dependence. No other assumptions were made on the environment and we expect such environments to display generic properties and to be representative of sufficiently general situations.
The model was studied analytically and numerically. The corresponding particle dynamics is given by an Itō differential equation with a multiplicative noise and no drift. We observed that random trajectories diffuse at large times. Also, the probability distribution of their displacements tends to a Gaussian at large times and the deviations from this asymptotic behavior decrease as t −1 . To take into account the fluctuations of mass due to the randomness of the environment, we introduced and studied the probability distribution of the density of mass. We obtained some analytical results on the density of mass in the case of a stationary environment and showed that it displays a power-law tail at large masses with exponent −2. In the general case, we observed competition between a trapping effect due to vanishing ejection rates of the random environment and the mixing due to its time dependence which leads to large fluctuations in the density of mass. These fluctuations were studied for both the smooth and the non-smooth random environment. In the smooth case, we showed that the PDF has an intermediate power-law behavior, like in the stationary case but with an exponent −3/2, followed by an exponential cutoff. Finally, we studied the spatial scaling properties of the mass distribution by introducing a coarse-grained density field. For smooth random environments the coarse-grained density was found to be scale invariant. We showed that at large masses, it possesses some scaling properties that depend on the coarse-graining size, the ejection rate typical amplitude and the Hölder exponent of the environment.
The overall dynamics of the model that was introduced in this paper contains several space scales and timescales. Mass rapidly accumulates near those regions with a vanishing ejection rate, and slowly moves following diffusion of the zeros. Depending on the properties of the environment, there is also a clear separation of length scales: small mass fluctuations at small scales and large ones at the scale of the distance between two vanishing ejection rates. This scale separation demonstrates the need to determine, by standard homogenization techniques, a large-scale effective diffusion tensor and possibly an effective transport term. We leave this topic for future work.
Extending our approach to dimensions higher than 1 is another possible direction for the future. The ejection rate will then vanish on complicated sets that, depending on the regularity, might display fractal properties. Varying the Hölder exponent of the random environment, its amplitude and the observation scales could then lead to a rather rich collection of different regimes. For instance, much more effort should be made to confirm whether or not a power-law is still present in the mass density probability distribution.
To conclude, we stress again that the presence of zeros in the environment has crucial effects on both the diffusive properties and the mass density statistics. They are responsible for mass accumulation but, at the same time, constitute barriers to transport. This duality implies that the fine statistical details of such systems crucially depend on the features of such zeros and in particular on the local structure of the ejection rate in their vicinity. In particular, the presence of an extra drift in the dynamics would completely alter the roles of the zeros and would prevent mass accumulation. In this work, we have decided to focus on ejection rates that can be written as the square of a generic Gaussian random function. However, this choice cannot always be relevant and, for instance, it might sometimes be more appropriate to write the ejection rate rather as the exponential of a random function. This would change drastically most of the results reported in this paper.