Passive advection of fractional Brownian motion by random layered flows

We study statistical properties of the process $Y(t)$ of a passive advection by quenched random layered flows in situations when the inter-layer transfer is governed by a fractional Brownian motion $X(t)$ with the Hurst index $H \in (0,1)$. We show that the disorder-averaged mean-squared displacement of the passive advection grows in the large time $t$ limit in proportion to $t^{2 - H}$, which defines a family of anomalous super-diffusions. We evaluate the disorder-averaged Wigner-Ville spectrum of the advection process $Y(t)$ and demonstrate that it has a rather unusual power-law form $1/f^{3 - H}$ with a characteristic exponent which exceed the value $2$. Our results also suggest that sample-to-sample fluctuations of the spectrum can be very important.


Introduction
The power spectral density of a time-dependent stochastic processes Y (t) is a meaningful feature of its spectral content which describes how its power is distributed over frequency. For stationary processes, it is usually defined as the time-average of the form where the symbol E {. . .} here and henceforth denotes the expected value with respect to different realizations of the process Y (t). For non-stationary processes, however, the expression given in (1) may not make sense and even the limit T → ∞ on the right hand side of it may not exist, so that one has to seek other meaningful interpretations of the power spectral density (see e.g., Refs. [1][2][3][4][5][6][7][8]). Although no unique tool exists for performing a time-dependent spectral analysis of non-stationary random functions, one of the physically plausible approaches [3,4] consists in using the time-averaged functional where W (f, t) is the Wigner-Ville spectrum The functions S f in (1) and W f in (2) and (3) become identical, once Y (t) is stationary. Many naturally occurring processes, as well as many processes encountered in engineering and technological sciences, exhibit power spectra of the form ∼ A/f α , where A is an f -independent amplitude. The characteristic exponent α may be as small as α = 1 in the case of flicker noise [9] or standard Sinai diffusion [10,11], α = 2 in the paradigmatic case of Brownian motion [6,9] ‡ and may even exceed the value of 2, e.g., for the Wigner-Ville power spectra in (2) and (3) of a super-diffusive fractional Brownian motion, for which one has α = 2H + 1 [4]. While the examples with α = 1 or α > 2 are rather rare, there exist numerous processes for which the power spectrum is a powerlaw function with 1 < α < 2. Few stray examples include electrical signals in vacuum tubes, semiconductor devices and metal films [1,9]. More generally, such a behaviour is observed in sequences of earthquakes [15] and weather data [16], in evolution [17], human cognition [18], network traffic [19], fractional Brownian motion with stochastic reset [20] and even in the distribution of loudness in musical recordings [21]. Recent experiments have also revealed such a behaviour of spectra for transport in individual ionic channels ‡ Note that while α = 2 is a valid result for Brownian motion, the observation of the A/f 2 law alone does not imply that one necessarily deals with the standard diffusion. The A/f 2 law holds for the finite-T power spectra in (1) of super-diffusive fractional Brownian motion with the Hurst index H > 1/2 [6], anomalous scaled diffusion [8], a variety of the so-called diffusing-diffusivity models [12] and also for the running maximum of Brownian motion [13] and diffusion in periodic Sinai chains [14], to name but a few processes. The amplitude A in the first three examples is, however, ageing, i.e., it is dependent on the observation time T . [22,23], electrochemical signals in nanoscale electrodes [24], bio-recognition processes [25] and intermittent quantum dots [5,26]. Many other examples and unresolved problems have been discussed in Refs. [5,[26][27][28][29].
In this paper we discuss a physical model with a super-diffusive dynamics and calculate its Wigner-Ville power spectral density defined in (2) and (3), which is shown here to exhibit a power-law with an exponent α exceeding 2. We focus on a stochastic passive advection process Y (t) in the presence of quenched, random layered flows. The model has been introduced originally by Dreizin and Dykhne [30] for the analysis of conductivity of inhomogeneous media in a strong magnetic field, and by Matheron and de Marsily [31] for the analysis of transport of solute in a stratified porous medium with flow parallel to the bedding. Subsequently, different facets of this model have been analysed in a great detail (see, e.g., Refs. [32][33][34][35][36][37][38][39][40]), establishing also a link, on a mathematical level, to the well-known models of statistical mechanics such as diffusion in presence of sources and sinks, spin depolarisation in random fields, self-repulsive polymers, and an electron in a random potential (see Ref. [36]). The model has been also generalised to study dynamics of more complicated objects, e.g., flexible polymers, subject to such velocity fields (see, e.g., Refs. [41][42][43][44][45][46]). While in the original settings in Refs. [30,31] a random motion X(t) in the direction perpendicular to the layered velocity fields was supposed to be diffusive, these latter works [41][42][43][44][45][46] provided some insight into the statistical properties of the process Y (t) for anomalous inter-layer diffusion. In particular, when X(t) is a standard diffusion, one finds the super-diffusive behaviour of the form E {Y 2 (t)} ∼ t 3/2 for the displacement along the direction of the flows, where the angle brackets here and henceforth denote averaging with respect to the distribution of the flow velocities. For a tagged bead of an infinite Rouse polymer, an even stronger super-diffusion has been predicted [41], E {Y 2 (t)} ∼ t 7/4 . Here, we consider a more general case when X(t) is a fractional Brownian motion with the Hurst exponent H, (0 < H < 1), and derive an H-parametrised family of super-diffusive laws for E {Y 2 (t)} and the corresponding family of the Wigner-Ville power spectra in (2), The paper is outlined as follows: In Sec. 2 we describe our model and introduce basic notations. In Sec. 3, we calculate the disorder-averaged mean-squared displacement along the Y -axis and also quantify its sample-to-sample fluctuations by analysing the coefficient γ v of variation of the corresponding distribution of E {Y 2 (t)}. We show that for the Hurst index H 0.22, the coefficient γ v of variation is less than unity, meaning that the sample-to-sample fluctuation are significant but their overall effect is not dramatic. On contrary, for 0 < H 0.22, γ v exceeds unity which implies that the standard deviation becomes larger than the mean value. Therefore, for such H, one expects the sample-to-sample variations of E {Y 2 (t)} to become much more important. Further on, in Sec. 4 we turn to the central point of our analysis -calculation of the disorder-averaged Wigner-Ville spectral density W f of the passive advection process. We show that the latter exhibits the frequency f -dependence of the form A/f 3−H , i.e., has an exponent α which exceeds the value 2. Here we also attempt to quantify the effective broadness of the distribution of the realisation-dependent value of W f . To this end, we consider the behaviour of W f for a finite observation time T and zero frequency, f = 0. We show that in this particular case, the coefficient of variation of the distribution grows (as a power law) with T , which signals that the sample-to-sample fluctuations are most important in this limit. Next, in Sec. 5 we concentrate on the representation of the Wigner-Ville spectrum in form of an Ising-type chain of "spins" σ k and analyse the behaviour of the effective time-dependent couplings J k,k ′ (t) between the "spins" σ k and σ k ′ . Finally, we conclude in Sec. 6 with a brief recapitulation of our results and outline the perspectives of further research.

The model
Consider the dynamics of a particle in the two-dimensional model system depicted in Fig.1, which we define as follows: a particle undergoes an unbiased random motion X(t), (starting at the origin at t = 0, X(0) = 0), along the X-axis. This random motion is the so-called fractional Brownian motion (fBm) [47]. The latter is formally defined as a stochastic integral with respect to the white noise measure dB s : where H ∈ (0, 1) is called the Hurst index. This is a Gaussian process with zero mean and covariance The standard Brownian motion with independent increments is recovered for H = 1/2. When H = 1/2, the increments are correlated so that for H > 1/2 if there is an increasing pattern in the previous steps, then it is likely that the current step will be increasing as well, resulting ultimately in a super-diffusive motion. For H < 1/2 the increments are negatively correlated, which entails a sub-diffusive motion. Further on, we mark along the X-axis the points X k = x k, k = 0, ±1, ±2, . . . (dashed lines in Fig. 1), where x is the distance between adjacent points. For simplicity, we set this distance equal to unity. At these points we have flows, depicted by arrows in Fig.1, with a velocity which is constant along the Y -axis (this constant is set equal to 1) and with k-dependent orientation described by a quenched random variable σ k = ±1, such that σ k = +1 corresponds to a flow in the positive Y -direction (positive velocity), and σ k = −1 corresponds to a flow orientated in the negative Y -direction. We concentrate here solely on the case when σ-s are delta-correlated and when there is no global flow in the system, i.e., when σ k ≡ 0.
Next, we assume that when the particle appears in-between the neighbouring points X k , it does not experience any drift in the Y -direction. However, once it arrives at any X k , it is passively (and instantaneously) advected on a fixed distance y = 1 along the Figure 1. A realisation of a pattern of random layered flows. The particle undergoes a fractional Brownian motion along the X-direction and is passively advected along the quenched random flows in the Y -direction.
direction defined by the arrow (flow) at this point. Consequently, the current position of the particle along the Y -axis obeys where N k [X(t)] is a trajectory X(t)-dependent random variable -the so-called local time [48,49] -which measures how many times within the time interval (0, t) the point X k has been visited by X(t). This random variable can be formally represented as where δ(. . .) denotes the Dirac delta-function.

Disorder-averaged mean-squared displacement of Y (t).
We analyse first the mean-squared displacement of the particle along the Y -axis. Thanks to (6), we can formally write this average as In order to calculate the expected value of the product of two local times taken at two different positions in space, we need to know the two-point probability distribution function for the fBm. The probability P (X(τ 1 )|τ 1 ; X(τ 2 )|τ 2 ; 0|0) that a fBm, starting at the origin, appears at position X(τ 1 ) at time moment τ 1 and at position X(τ 2 ) at time moment τ 2 (τ 1 and τ 2 being unordered) is given by where g is the correlation coefficient Since Therefore, by using (9) we readily find that The above expression represents the desired two-point correlation function of the local times of a fBm at two different (or coinciding) points taken at the same time moment t. By inserting (11) into (8), performing the averaging over the distribution of {σ k } and summation over k, (as well as appropriately changing the integration variables), we get where and θ 3 (. . .) is the Jacobi's theta function. Turning to the limit t → ∞, we find where the symbol o t H means that the omitted terms grow with t slower than t H . Note that taking into account only the leading in the limit t → ∞ term in the right hand side of (14), is tantamount to converting the summation over k into an integral -by using the Euler-Maclaurin summation formula -and discarding all the correction terms. Consequently, we arrive at the following asymptotic large-t form which represents the desired result on the H-parametrised family of anomalous laws describing the passive advection of a fBm by random layered flows. Note that for any H ∈ (0, 1) the disorder-averaged mean-square displacement exhibits a super-diffusive behaviour. In particular, for H = 1/2 we recover the result E {Y 2 (t)} ∼ t 3/2 obtained previously in Refs. [30,31]. For H = 1/4, we find an even stronger super-diffusive law E {Y 2 (t)} ∼ t 7/4 , obtained earlier in Ref. [41] for the dynamics of a tagged monomer in an infinite Rouse polymer in the presence of such random flows. In general, we observe that the smaller H is, the more pronounced becomes the super-diffusive behaviour of Y (t), which is a simple consequence of the fact that for a more spatially confined subdiffusion, the particle spends more time in a given layer, and hence, is advected on larger scales.
3.1. Sample-to-sample fluctuations of the mean-squared displacement.
Here we address the question of the sample-to-sample fluctuations of E {Y 2 (t)} as defined in (8) for different realisations of the patterns of random flows. In order to quantify these fluctuations, we focus on the variance of E {Y 2 (t)}, which is defined as follows We turn to the limit t → ∞ and concentrate on the leading term which dominates the large-time asymptotic behaviour of the variance. It is plausible then to convert the infinite sums in (16) into integrals, which can be readily made dimensionless by an appropriate change of the integration variables. In doing so, we get with γ v being the (time-independent) coefficient of variation of the distribution A straightforward calculation allows us to write the following expression for γ v : where q, for a given H, is a numerical factor defined by The coefficient of variation γ v is a meaningful characteristic of the effective broadness of the distribution showing how E {Y 2 (t)} fluctuates from sample to sample. Being unable to perform exactly the four-fold integral (19) entering the expression for γ v , we evaluate it numerically and plot it as a function of the Hurst index in Fig. 2. We observe that γ v is a monotonically decreasing function of H which perfectly makes sense because for a fixed time interval (0, t) the span of the trajectory X(t) is a monotonically increasing function of H. Further on, we notice that γ v (H) is less than unity for most of the values of H in the interval (0, 1) meaning that sample-to-sample fluctuations here are not very significant. The coefficient of variation exceeds slightly 1 (such that the standard deviation becomes greater than the mean value) only for H 0.22; note that γ v (0) = (4/3) 1/4 ≈ 1.075. Only in this region, i.e., for very spatially confined sub-diffusive processes X(t), the sample-to-sample fluctuations may become rather significant (see, e.g., Ref. [50,51]).
We turn next to the calculation of the corresponding power spectral density of the passive advection process Y (t). To this end, it is first expedient to somewhat simplify these expressions. We formally rewrite the covariance of the process Y (t) as with The Wigner-Ville spectrum of the process Y (t) then reads where θ(. . .) denotes the Heaviside theta-function. In (23), we can straightforwardly perform the integral over τ to get Inserting the latter expression into (3) and performing the integration over t, we find The expression in the right-hand-side of (25) is still quite complicated since it involves oscillating functions of τ 1 and τ 2 and the integration limits are mixed. We notice then that the existence of the limit in the right-hand-side of (25) implies that the two-fold integralˆ∞ grows linearly with T when T → ∞. In turn, this means that the Laplace transform of the expression in (26), (with respect to T with the Laplace parameter λ > 0), diverges as 1/λ 2 in the limit λ → 0. This permits us to formally rewrite (25) as which appears to be more convenient for further calculations. We focus next on the disorder-averaged Wigner-Ville spectrum W f . The disorderaveraged function Q(τ 1 , τ 2 ) reads Since we are interested the behaviour in the limit λ → 0, we have to focus on the asymptotic behaviour of the Jacobi theta-function in the limit τ 1 → ∞ and τ 2 → ∞.
In this limit the expression in the exponential in the theta-function tends to zero such that we find that, in the leading order, the theta-function behaves as which yields Inserting this expression into the two-fold integral in (27) and performing the integrations, we find where the symbol O (1) means that the omitted terms are constant in the limit λ → 0. Inserting the latter expression into (27) and taking subsequently the limit λ → 0, we obtain the following result for the disorder-averaged Wigner-Ville power spectral density: This result is the central point of our analysis and defines an H-parametrised family of super-diffusive spectra characterised by the exponent α = 3 − H > 2 for any 0 < H < 1. Notice also that the H-dependent overall factor is a positive definite, monotonically increasing function of H. We further note that for the case when X(t) is a standard Brownian motion, we have while for H = 1/4 we get which particular case corresponds to the dynamics of a tagged bead in an infinite Rouse chain.

Sample-to-sample fluctuations of the Wigner-Ville spectrum at zero-frequency
Consider the following finite-time averaged Wigner-Ville spectrum in (2) and (3) at a finite observation time T (such that we drop the limit T → ∞ in (2)) and denote the resulting expression as W (T ) f . Then, the disorder-averaged W with Q(τ 1 , τ 2 ) given by (30). In the limit f → 0 the integrand in the latter expression stays finite. Performing next the change of the integration variables τ 1 = T x 1 , τ 2 = T x 2 , we realise that the asymptotic behaviour of the finite-time disorder-average Wigner-Ville spectrum follows where Ψ 1 (H) is the function .
The above summations can be carried out precisely as it was already done in our previous calculations which lead us to (28). Consider, for the sake of simplicity, the first sum on the right hand side of (41). It is straightforward to show that where and Other sums can be tackled in essentially the same way. Noticing next that in the limit f → 0 the integrand in (40) stays finite, we find Similarly to the analysis of the mean value (36), here we perform the following change of the integration variables: τ j = T x j , j = 1, . . . , 4, which permits us to cast the variance into the form with Ψ 2 (H) being a dimensionless function of the Hurst index H. Combining (36) and (46), we conclude that the probability density characterising the finite-time Wigner-Ville spectrum at zero-frequency has a coefficient of variation γ W ∼ T 2−H/2 , implying that the distribution broadens with T and the sample-to-sample fluctuations become substantially more important. However, we are not in position to determine the coefficient of variation for finite f > 0 here. In principle, we expect that similarly to what happens with a super-diffusive fBm, γ v will attain a finite value in the limit T → ∞ (see [7]). This analysis goes beyond the scope of the present work and will be published elsewhere.

Effective couplings in the Ising-like model representation of the Wigner-Ville spectrum
In this last section we represent the Wigner-Ville spectrum W f as a Hamiltonian of an Ising-like chain of "spin" variables σ k , (which prescribe the directions of the flows), and analyse the form of the effective couplings in this representation. Taking advantage of our previous results, we have that the Laplace-transformed over the observation time T power spectral density admits the following form where the couplings J k,k ′ (λ) are defined by By noticing that J k,k ′ (λ) are evidently real-valued, even functions of f , we can also rewrite (47) in the form which is somewhat easier to handle. Still, we are able to determine J k,k ′ (λ) for arbitrary k and k ′ only for the Brownian motion with H = 1/2. For arbitrary H, it turns possible only to determine J k,k (λ).

Particular case H = 1/2
It is expedient to get first some idea on the form of J k,k ′ (λ). This can be done in the special case H = 1/2 in which the integrals in (48) can be performed exactly. After straightforward calculations, we find that for H = 1/2 the couplings in (46) are given explicitly by which expression holds for any value of the parameters f , λ, k and k ′ . Note that J k,k ′ (λ) are real-valued functions and we have chosen the form in (49), which involves the unit imaginary number i just for the sake of compactness. Note, as well, that J k,k ′ (λ) are oscillatory functions of |k − k ′ | with a period of oscillations ∼ 1/ √ f in the limit λ → 0. In Fig. 3 we show the couplings J k,k ′ (λ) as functions of k and k ′ for a fixed frequency f . Quite generally, J k,k ′ (λ) reach their maximal values at the origin, i.e., for k = k ′ = 0, and vanish for large k and k ′ . It is evident from Fig 3 that the decay of the couplings is spatially anisotropic, with k = k ′ being the direction of a slow decay and k = −k ′ -the direction of a fast decay. Upon decreasing the frequency, the spatial dependence of the couplings acquires further anisotropic structure, as illustrated in the right panel of  Several points have to be emphasised. First, we notice that the result in (49) permits us to reproduce (32) in the particular case H = 1/2. Indeed, noticing that σ k σ k ′ = δ k,k ′ , we have which is precisely our result in (33), obtained by taking into account only the leading asymptotic behaviour of the Jacobi theta-function. This means, in turn, that the arguments used for the derivation of the result in (33) are correct.
Second, we notice that taking the limit λ → 0 and performing the summation operations can not be interchanged. Indeed, for any fixed k and k ′ we have lim λ→0 J k,k ′ (λ) = 0, which signifies, in turn, that for small λ the sums in (46) are dominated by the terms with large k and k ′ . Further, since J k,k ′ (λ) are linked to the correlation function of the local occupation times at two different points, they depend simultaneously on both the distances |k| and |k ′ | from this sites to the origin (the starting point of the trajectory X(t)) and also on their relative distance |k −k ′ |. In these dependences, which are simple exponential functions for H = 1/2, the characteristic decay of J k,k ′ (λ) with |k − k ′ | depends on both f and λ and here we can safely set λ = 0. On contrary, we have to keep the dependence of the characteristic decay length on λ in the terms dependent only on k and k ′ (the sum of two exponential functions in the first factor in the right hand-side of (49).
We turn next to the general case of an arbitrary H ∈ (0, 1) focussing on the diagonal terms. For k = k ′ , the couplings J k,k (λ) in (48) read which form often arises in the analysis of various physical problems via the optimal fluctuation method (see, e.g., Ref. [52,53]). On the other hand, the integral in (55) defines the so-called special Krätzel function (see, e.g., Ref. [54]) and therefore, the diagonal couplings can be formally represented as Before we proceed to the analysis of the asymptotic behaviour of J k,k (λ) in (57), it seems expedient to consider first a particular case of Brownian motion (H = 1/2). Here, we have such that which expression agrees perfectly well with our previous result in (49) with k set equal to k ′ . Further on, summing the expression in the second line in (55) over all k, and taking an appropriate limit in the resulting Jacobi theta-function, we recover exactly our result in (32) for the disorder-averaged Wigner-Ville spectrum.
The k-dependence of the diagonal couplings given in (57) is plotted in Fig. 4. The case of the Brownian motion, corresponding to (49) with k = k ′ , is depicted with a dotted line. We observe that upon increasing the Hurst index, the diagonal coupling increases monotonically for any k, while for any fixed H the large-|k| behavior indicates the stretched-exponential decay predicted by the asymptotic result (63).
Lastly, we consider the asymptotic behaviour of J k,k in (57) in the limits ≡ |k|λ H ≪ 1 and |k|λ H → ∞. For the sake of notational convenience, we rescale k by defining κ = kλ H / √ 2.

Conclusions
To recap, we have studied here a model of anomalous diffusion of a tracer particle in presence of stratified layers of quenched random flows. The particle undergoes a fractional Brownian (fBm) motion with the Hurst index H ∈ (0, 1) in the direction perpendicular to the flows (along the X-axis) and is passively advected in random directions (with no global bias) by the flows along the Y -axis. Averaging the squared displacement Y 2 (t) along the flows over the flow realisations, as well as over the fBmtrajectories, we obtained the mean-squared displacement along the flow direction and showed that for large times it grows in proportion to t 2−H , which law defines a family of super-diffusive processes. In order to quantify the sample-to-sample fluctuations of the displacement along the random flow, we computed the coefficient of variation of the probability distribution function of the displacement. We showed that such fluctuations are typically not very important for H 0.22 but may become important for sufficiently small values of H. Next, we introduced the Wigner-Ville power spectral density (see (2) and (3)) of the random process Y (t) and derived an exact result showing that the disorder-averaged Wigner-Ville spectrum < W f > scales with the frequency as f H−3 , which defines a family of rather unusual super-diffusive spectra. To estimate the sample-to-sample fluctuations of the spectrum, we examined the limit of zero frequency of the Wigner-Ville spectrum averaged over a finite time interval of duration T . We found that the ratio of the standard deviation of the spectrum and of the mean value diverges with T , which implies that the distribution of the spectrum with respect to different realisations of disorder progressively broadens with the observation time.
Lastly, representing the Wigner-Ville spectrum of the process Y (t) for a given realisation of flows as an Ising-like model of "spins" σ k , which define the local direction of the flows along the X-axis, we have analytically determined the coupling terms in this model, in the Laplace domain. In particular, we showed that the latter exhibit a non-trivial dependence on the distance from the origin and the Laplace parameter.
In general, our analysis reveals that sample-to-sample fluctuations are important, both for the dynamical characteristics and the spectral behaviour. This issue is beyond the scope of the current work and will be examined elsewhere.