Universal and Non-Universal Features in the Random Shear Model

The stochastic transport of particles in a disordered two-dimensional layered medium, driven by correlated y-dependent random velocity fields is usually referred to as random shear model. This model exhibits a superdiffusive behavior in the x direction ascribable to the statistical properties of the disorder advection field. By introducing layered random amplitude with a power-law discrete spectrum, the analytical expressions for the space and time velocity correlation functions, together with those of the position moments, are derived by means of two distinct averaging procedures. In the case of quenched disorder, the average is performed over an ensemble of uniformly spaced initial conditions: albeit the strong sample-to-sample fluctuations, and universality appears in the time scaling of the even moments. Such universality is exhibited in the scaling of the moments averaged over the disorder configurations. The non-universal scaling form of the no-disorder symmetric or asymmetric advection fields is also derived.


Introduction
In 1973, Dreizin and Dykhne proposed a model for the density of a diffusion material in stationary current flow, through a disordered conducting medium [1]. In particular, the model was conceived for calculating the effective conductivity of a nonuniform concentration of carriers, describing its stationary transport under the influence of both diffusion and convection (with velocity U). The carriers transport was assumed to be two-dimensional and limited to a stripe of transverse dimension L, see Figure 1. As the fluctuations of the symmetric conductivity tensor in the transverse direction were taken to be large compared with longitudinal components, the transverse diffusion overwhelms the longitudinal one, i.e., D ≡ D y D x . In the transverse (x) direction, the diffusing particles experience random pulsations, which produce narrow convective flow with a velocity U. Here, the velocity field has a "quenched" randomness, from which the expression random shear model arises. Therefore, the stochastic motion of a generic tracer inside the channel is a combination of two mechanisms. On one side, the tracer undergoes normal diffusion on y, the channel's transverse direction, on the other, it experiences a random static velocity field directed toward the channel's axis x. In terms of Langevin equations for an ensemble of independent particles, the model is hereby defined aṡ where{ξ i (t)} are independent, delta-correlated and zero-mean Gaussian thermal noises. Periodic boundary conditions are enforced on the y-direction at y = ±L/2, to implement a channel-like geometry along the x-axis. The above equations were introduced in 1980 by Matheron and de Marsily [2] in their celebrated paper on the solutes transport in porous media. Since its introduction, the most striking feature of the random shear model was its longitudinal superdiffusive behavior, δ 2 x(t) ∼ t 3/2 . Generally speaking, the average · · · is a twofold operation, namely, an average over the thermal noise in (2) and an average over the ensemble of the disorder layered velocity field. In the late-80s-early-90s period, there was a flourish of mathematical models devoted to the investigation of the Drezin-Dykhne superdiffusive regime and its relation with statistical properties of the quenched advection field. The vast majority of the theoretical and numerical studies assumed the distribution P(U) of the static shear field to be Gaussian, with a stationary autocorrelation function in space, of the type U(y 1 )U(y 2 ) ∼ δ(y 1 − y 2 ) [3][4][5][6][7][8][9][10][11][12]. However, although this simplified picture was largely adopted, the necessary and sufficient conditions for having a superdiffusion ∝ t 3/2 were already clearly stated in the Dreizin-Dykhne paper: i) the stationarity of the advection field autocorrelation function U(y 1 )U(y 2 ) = f (|y 1 − y 2 |), and ii) its integrability ∞ −∞ dy U(0)U(y) < ∞ [1]. Indeed, Matheron and de Marsily assumed a stationary Gaussian-shaped velocity space autocorrelation, a property which was later used in the analysis reported in References [13][14][15]. However, models releasing either one of the aforementioned hypotheses were also considered. As an example, a stationary although power-law velocity autocorrelation function was proposed in Reference [16] to account for solutes diffusion in certain types of fractured rocks, while the absence of stationarity characterizes the velocity autocorrelation function in Reference [17]. As anticipated by the analysis of Dreizin and Dykhne, any deviation from the stationary short-range velocity autocorrelation function entails a superdiffusive behavior different from the classical t 3/2 law. In general, the statistical properties of the disordered velocity field are limited to the Eulerian two-point autocorrelation function, with the only exception of the work in [18] which introduced an n-point correlation function of the type U(y 1 )U(y 2 ) · · · U(y n ) = σ n f n (y 1 − y 2 ) f n (y 3 − y 4 ) · · · f n (y n−1 − y n ), and also substituted the Brownian dynamics along y, Equation (2), with the anomalous fractional diffusion. Furthermore, shear models with a deterministic dependence of the velocity on the coordinate y were also considered, together with their impact on the time behavior of the mean square displacement: examples are a linear advecting field U(y) = U 0 y [19][20][21], a power-law U(y) = U 0 |y| β sign(y) [22], or a sinusoidal U(y) = U 0 sin(y) [23].
The entire body of the theoretical works devoted to the random shear model takes into consideration the average of any observable on both the thermal noise and the disorder advection velocity field. Usually, the scaling of the solutes probability density function (PDF) along the channel's axis x, P(x, t), and of its moments, [x(t) − x(0] m , are studied within this framework [3,5,11]. The interest stems from the typical experimental condition where a stationary flux of water is flowing along a preferential direction, say x, with random static orientation, across an ensemble of disordered media, such as several porous samples or fractured rocks, see Figure 1. Imagine now dropping a tracer inside one of these porous samples and observing its dispersion in time. Repeating this transport experiment several times on the same sample mimics the thermal average of the tracer's diffusion. Performing an additive average over the porous media in the ensemble gives the double average considered in previous works. In Reference [24], we adopted another point of view, which was also partially embraced in Reference [15]. We have calculated the tracers' PDF in the case of a particular case of symmetric quenched disorder, performing the average over a set of uniformly spaced initial conditions along y. This point of view originates from the experimental perspective of having at hand one only sample, namely a single fractured rock or, in general, a single disordered medium, on which the diffusion experiment can be conducted several times. However, in this case a large amount of independent tracers is released inside the channel, a solute density uniformly spread throughout the channel breadth, but at the same position x. The diffusion experiment consists in observing how the initial density is dispersed in x and in time by the water flux. This corresponds to performing an average over the initial conditions. Repeating this experiment several times yields the statistical or thermal average. In this paper, we carry out a more thorough and extensive study on both the velocity field statistical properties and the position moments in the random shear model. In doing so, we unveil how the scaling properties exhibited by the velocity field are intimately connected to the spreading in time of the position moments. Generally speaking, the computation of the moments is used to test the scaling form of the propagator P(x, t) in the large-time limit [5,11,24]. In line with the model proposed in Reference [24], we perform our theoretical and numerical analysis by introducing a long-ranged velocity field whose properties can be fully established. We derive a simple and precise scaling framework for moments averaged over the disorder, and for moments averaged over an ensemble of uniformly spaced initial conditions, given a quenched configuration of the advection velocity field. We demonstrate the presence of universal features and lack of universality due to large sample-to-sample fluctuations. Importantly, we analyze and discuss the scaling of the moments in the cases of zero disorder, i.e., for deterministic shear models with long-ranged velocity fields. We show that all the scattered results present in the literature on this subject fit neatly into our mathematical scheme.

Velocity Field Properties
We discuss the static and dynamical properties of a shear longitudinal field which, according to References [23,24], is generated by a superposition of M sinusoidal modes, or cosinusoidal modes: |k n | γ/2 cos(k n y + φ n ).
The wave numbers k n run over the set k n = 2π/L(1, 2, . . . , M) and −1 < γ < 1. The Equations (3) and (4) can be expressed in a compact form as with the following prescriptions: (3) is valid, and for Equation (4), with U k 0 = 0. This choice carries three important advantages when compared to the previous definitions of the velocity field implemented in former works. First, it is automatically verified that Second, the way the disorder is implemented is through the phases φ n (n = 1 . . . , M), each one drawn from a uniform independent distribution on the support [0, 2π]. Third, all the statistical properties of the advection shear field can be easily calculated. For instance, the correlation function appears to be which yields a discrete spectrum ∼ k γ n . Therefore our model recovers the cases studied in Reference [16], where the shear field was introduced through its first (zero) and second (power-law) moments. However, in our formulation we overcome any divergence problem, thanks to the regularization due to the sum over discrete k n . Moreover if U k n =const, in the limit of M → ∞, our formulation recovers the Dikhne model [3][4][5][6][7][8][9][10][11][12].

Static Properties
The four points static correlation function is which is different from the expression proposed in [18]. Importantly, Equation (10) implies U(y) 4 φ = 3 U(y) 2 2 φ , and it can be shown that any even moment enjoys the Gaussian property U(y) m φ = (m − 1)!! U(y) 2 m/2 φ , while the odd moments vanish. This means that if we look at the probability of having a shear U at a given channel's position y, over several disorder configurations {φ n }, i.e., P φ (U|y), this would be in the average a Gaussian. In addition, since U(y) is a stationary process in space, its PDF must be the same for any position y. This is indeed shown in Figure   The collected values of U(y) involve 2000 uniformly distributed points along the y axis. The red curve represents the same PDF for the velocity field U(y) sampled over 10 5 points along the channel's width. The other simulation parameters coincide with those of panel (a). Panel (c): same quantity as in panel (b) obtained from the velocity field in the inset, namely with U(y) given by (3) and φ n = 0, ∀n. The shear field in the inset shows the presence of antisymmetric long jets close to y = 0 and, smaller values elsewhere. Other simulation parameters coincide with those of panel (a,b). The red curve represent the Gaussian P(U) shown in panels (a,b). Panel (d): same quantity as in panel (b,c) obtained from the velocity field in the inset, namely with U(y) given by (4) and φ n = 0, ∀n. In the inset the shear field is dominated by symmetric shoots close to the origin (U(−y) = U(y)) surrounded by negative values of the velocity U(y). The red curve represent the Gaussian P(U) shown in panels (a,b). Other simulation parameters coincide with those of panel (a-c).
Let us now take a different point of view that, in our opinion, has been overlooked in previous works. Let us consider the velocity field when the disorder is quenched: this corresponds to fixing the values of the phases {φ n } in (5) 2πσ 2 and this is true for almost any random configuration of the disorder, i.e., any typical {φ n }.
However, it is also true that for specific choices of the disorder, the distribution P y (U|φ) may markedly differ from a Gaussian. In Reference [24], we addressed this problem by considering a shear field of the type (3) and quenched phase configuration {φ n }, where the single phases φ n could be either 0 or π with probability 1/2. In this paper, we consider a different case, namely the complete absence of disorder, i.e., φ n = 0 ∀n. Specifically, we will consider two situations: an ordered shear field (3) (the case M = 1 has been proposed in [23], as already mentioned), and an ordered shear field (4). The y-dependence of these two fields are reported in the insets of panels (c) and (d) of Figure 2, respectively. It is possible to see that, although the asymmetry of the shear field (3) U(−y) = −U(y) guarantees the symmetry of P y (U|φ = 0), on the contrary, the parity of the cosine (U(−y) = U(y)) entails an asymmetric PDF, a situation which was already mentioned in Reference [22] within the context of deterministic power-law shear. Importantly, both the shear field distributions, P y (U|φ = 0) are not Gaussian, as it can be appreciated in the main panels (c) and (d) of Figure 2.

Dynamical Properties
From the analysis reported above, it appears evident that two types of averages can be performed independently: one is taking a well-defined height y on the channel and studying any observable by averaging over the disorder, i.e., over the ensemble of the quenched random phases {φ n }. This is the point of view usually taken in previous works [1,2,5,6,[8][9][10][11][12]14,16,22,25,26]. The other point of view is to quench the disorder and subsequently average the observable over its dependence on y ∈ [−L/2, L/2]. This kind of analysis has been adopted in [24] and partially in [15], and it needs to be handled carefully when taking into account dynamical observables. As a matter of fact, the solution of the stochastic Equation (2) is where w t indicates a Wiener's process, i.e., w t = 0 and w s w t = |t − s|. Hence, the time dependence of the shear flow (5) is expressed as implying that, for a dynamical observable, the average over y coincides with averaging over a uniform set of initial conditions y 0 , i.e., For instance, from (11) the mean drift of a particle is zero whether the average is taken over the initial conditions or over the disorder, i.e., U(t) w φ = U(t) w 0 = 0, where . . . w represents the average over the thermal noise. While analyzing the two-time correlation function, we first perform the thermal average: This expression clarifies how, without any further averaging procedure, the velocity process is neither stationary, nor satisfies the independence on the initial condition and, on top of that, it shows an explicit dependence on the particular choice of the quenched disorder. Now, performing an average over the phases, one easily gets for the time-ordered situation in which t 1 > t 2 . Thanks to the fact that e i(k n 1 +k n 2 )y 0 0 = e i(φ n 1 +φ n 2 ) φ = δ(n 1 − n 2 ), the same result is obtained performing the average over the initial conditions, i.e., U(t 1 )U(t 2 ) w φ = U(t 1 )U(t 2 ) w 0 . The equation (12) is stationary and independent both from the initial conditions and from the specific choice of the disorder taken into account. Moreover, Equation (12) can also be achieved by applying the following formula [5,[9][10][11] where P(y, t) represents the one-dimensional Brownian propagator on a strip of amplitude L with periodic boundary conditions, i.e., Let us now turn to the study of the three (ordered) time velocity correlation functions. The disorder averaged correlation function U(t 1 )U(t 2 )U(t 3 ) w φ is easily shown to be 0. On the contrary, the average over the initial positions yields From the properties (6) and (7), it appears clear that taking an asymmetric velocity field such as (3), or its symmetric counterpart (4), yields in principle different results. After straightforward passages the previous expression becomes This correlation function exhibits an apparent dependence on the particular choice of the quenched disorder, which makes this observable unpredictably oscillating around 0. Moreover, since it depends only on the time difference t i − t j (t i > t j ), it highlights that the shear field (averaged over the initial conditions) is a strict sense stationary (SSS) process [27]. Thanks to this property, the expression (17), derived under the condition t 1 > t 2 > t 3 , is equivalent to those obtained for any other time ordering by a simple permutation of the indices [27,28].
The disorder-averaged four-time correlation function, when t 1 > t 2 > t 3 > t 4 , is defined as 2 4 e i(φ n 1 +φ n 2 +φ n 3 +φ n 4 ) φ e i √ 2D(k n 1 w t 1 +k n 2 w t 2 +k n 3 w t 3 +k n 4 w t 4 ) w e i(k n 1 +k n 2 +k n 3 +k n 4 )y 0 , which is equivalent in both cases (3) and (4). Thanks to the identity e i(φ n 1 +φ n 2 +φ n 3 +φ n 4 ) φ = δ(n 1 + n 2 )δ(n 3 + n 4 ) + δ(n 1 + n 3 )δ(n 2 + n 4 ) + δ(n 2 + n 3 )δ(n 1 + n 4 ), the final expression reads As in the case of Equation (14), the same result can be obtained by applying the formula Expression (19) allows drawing the conclusion that the average over the phases entails the disappearance of the y 0 -dependence, a characteristic which is also evident in Equation (13). Moreover, the explicit dependence on the time difference t i − t j (t i > t j ) also makes the disorder-averaged velocity a SSS. Hence, U(t 1 )U(t 2 )U(t 3 )U(t 4 ) w φ is invariant to the permutations of t 1 , t 2 , t 3 and t 4 . We now study the four-time correlation function averaged over the initial conditions, i.e., 2 4 e i(k n 1 +k n 2 +k n 3 +k n 4 )y 0 0 e i √ 2D(k n 1 w t 1 +k n 2 w t 2 +k n 3 w t 3 +k n 4 w t 4 ) w e i(φ n 1 +φ n 2 +φ n 3 +φ n 4 ) .
A direct calculation yields .
Once again, like in Equation (17), the explicit dependence on the quenched phase configurations is exhibited.
In view of the above analysis, we can draw the following conclusions about the shear field statistical properties:

•
The shear field is a SSS, both averaged over the disorder and over a uniform distribution of initial conditions

Position Moments
Owing to the properties investigated in the former section, in the present we show the explicit calculation of the second, third and fourth moment averaged over the phases [x(t) − x(0)] m w φ and over the initial conditions [x(t) − x(0)] m ] w 0 . The general formulas that we will adopt are the following [9][10][11]25]

Second Moment
While the drift (m = 1) is invariably 0, the second moment is the only moment that has the same expression if averaged over the disorder or over a uniform distribution of initial conditions. This is true for any configuration of the phases, including the no-disorder cases, as well as for both choices (3) and (4). It can be easily deduced from the definitions (23) and (24) by using the correlation function (13). Hence the time behavior of can be studied through the Laplace transform The pole at 0 dictates the long time behavior ∼ Lt, while for short and intermediate times t The correct analysis can be performed by inverting (25) [ and is reported in Reference [24]. We hereby present only the scaling behaviors (see Figure 3): One sees that the asymptotic normal behavior is restored thanks to the presence of the channel's boundaries, which impose an exponential cut-off to the velocity autocorrelation function (13) [24]. As anticipated, the unique feature of the mean square displacement (26) is that [x(t) − x(0)] 2 w 0 has the same expression for whatsoever disorder quenched configuration, as it is shown in Figure 3 (blue symbols), and this corresponds exactly to the disorder-averaged MSD [x(t) − x(0)] 2 w φ (solid lines).

Third Moment
The position third moment is identically zero in the case of disorder average. On the other hand, it is non-zero if the average over the initial conditions is performed, i.e., where the velocity correlation function is given by (17). Tracing the same procedure of the second moment, we first calculate its Laplace transform At first glance, one sees that the third moment carries an explicit dependence on the disorder, in analogy to (17). Hence, since a generic disorder does not preserve the symmetry x → −x, the third moment cannot be neglected. As a matter of fact, in Figure 4a   A well-defined non-zero scaling in time appears manifest only by considering the symmetric shear case with zero disorder, namely the field (4) with {φ n } = 0, ∀n. A rather lengthy calculation, briefly reported in the Supplementary Information, shows that the three scaling regimes of the third moment can be identified on the score of what we have conducted for the second moment. We hereby report the overall result

Fourth Moment
We analyze the case of the fourth moment, starting with the disorder-averaged case where the expression (19) must be inserted. Performing the Laplace transform of the former equation consists of transforming each of the terms appearing in (19). In the Supplementary Information, we detail the procedure to highlight the time scaling expression of [x(t) − x(0)] 4 w φ , which retraces that of the second and third moments reported in the previous sections. We sum up the result of our analysis in the following formula The analysis of the fourth moment averaged over the initial conditions is more involved as it stems from the formula A close inspection of the velocity correlation function (22) reveals that when the conditions n 1 = −n 2 , n 1 = −n 3 and n 2 = −n 3 are satisfied, Therefore, we can split the fourth moment, averaged over the initial conditions, into two overall contributions: ∑ n 1 ,n 2 ,n 3 =−M n 1 =−n 2 n 1 =−n 3 n 2 =−n 3 f (φ n 1 , φ n 2 , φ n 3 , φ n 1 −n 2 −n 3 ; t).
The last term represents the explicit dependence of [x(t) − x(0)] 4 w 0 upon the quenched configuration of the disorder. For a generic random choice of {φ n } and large M we expect that this term does not give a substantial contribution to the scaling behavior of (32). Indeed, in Figure 4b, several curves pertaining to random configurations of {φ n } exhibit small fluctuations around the average value represented by [ (31)) for almost any random typical disorder.
On the other hand, if the disorder is quenched to zero, namely if the velocity advection field is (3) or (4) with {φ n } = 0 ∀n, the second term in the Equation (32) plays a significant role, leading to scaling forms different than those reported in (31). These are Supplementary Information reports the calculations needed to derive the expression (33).

General Scaling of the Moments
We now generalize the results obtained in the previous section, designing a complete scaling scheme for the moments of order m in the three time regimes, when different advection field conditions are analyzed.
Generally speaking, when the phase disorder is quenched, and for M → ∞, the overall majority of the systems will display apparent scaling behavior only for even moments, while the odd moments show large sample-to-sample fluctuations around zero. Thus, for a generic quenched disorder, the even moments averaged over the initial conditions display a typical scaling, i.e., compatible with the scaling of moments averaged over a statistical ensemble of the phases configurations (average over the disorder). The cases where the disorder is absent, i.e., when φ n = 0 ∀n, are not typical. They are indeed characterized by moment scalings in time different from those obtained from the moments averaged over the disorder. Moreover, while the asymmetric case of sine advection field (3) does not possess well-defined odd moments scaling, this can be evaluated explicitly in the case of symmetric shear velocity field (4).

Disorder-Average and Random Phase Configuration
A direct, although tedious, calculation shows that the scaling of the moments averaged over the disorder is Here, the index m can assume only even values, as the odd moments are identically zero. The expression (34) is confirmed by the numerical simulations in Figure 5  1/m φ yields a perfect collapse of the different moments in each of the three regimes. This is a bright indication that the even moments exhibit an anomalous scaling [29][30][31] in the three regimes.
We have seen in the Sections 3.2 and 3.3 ( Figure 4) that the third and fourth moments, averaged over the initial conditions, exhibit the average trend of [x(t) − x(0)] m ] w φ . However, for a generic choice of the phases, fluctuations around this mean appear to increase with time. While these sample-to-sample fluctuations are considerably marked in case of odd moments, wildly scattered around zero, when m is an even number they do not appreciably affect the scaling in (34). This is indeed shown in Figure 5, where a single realization of [x(t) − x(0)] m ] w 0 obtained for a random choice of {φ n } (colored symbols) overlaps with the disorder average moments (34).

No-Disorder: Symmetric and Asymmetric Quenched Velocity Fields
When the no-disorder condition is enforced, i.e., φ n = 0 ∀n, the scaling of the moments differs substantially from the disorder-averaged case (34). A lengthy calculation, although conceptually non-demanding and equivalent to those reported in Supplementary Information, shows that, when the symmetric velocity advection field (4) is chosen, even and odd moments exhibit clear scaling forms in time (see solid lines in Figure 6). In particular, it appears that these scaling expressions have a continuous linear dependence on m in the range of short to intermediate times t  When the velocity field is asymmetric, i.e., (3), the odd moments present large fluctuations in time, although they are zero on average. Therefore, a precise scaling form in time cannot be inferred. The even moments instead follow the same scaling form of the symmetric case in the short to intermediate time regime t L 2 D , and they coincide for large times (dashed lines in Figure 6).
Summing up, the results for the symmetric and asymmetric advection fields are hereby reported The former Equation (35) holds for any m when the velocity field is symmetric, namely in the cosine case (4). Its validity is limited to only values of 2m for antisymmetric velocity fields (3).

Conclusions
We have furnished a thorough analysis of the random shear model in the case of long-ranged advection fields. The strength of the correlations along y is dictated by the exponent of the discrete spectrum of the velocity field ∼ k γ n . The elegant way in which the random shear velocity is introduced presents several advantages if compared to the previous analysis, among which we can enumerate the possibility to derive exactly all the advection (shear) field statistical properties by means of the n−points correlation functions and of the PDF P(U). Importantly, it allows the straightforward implementation of two kinds of averaging procedures, one · · · w φ over the thermal noise w and the shear field disorder embodied by the phases {φ n }, the other · · · w 0 over the thermal noise w and a uniform set of initial conditions y 0 spaced along the channel's size L, for a quenched configuration of the disorder. In particular, we focused on the comparison between the scaling form in time of the position moments achieved by averaging over the disorder ensemble, and the scaling exhibited by the moments averaged over the initial conditions. Three distinct cases of the quenched advection field were thoroughly analyzed: (i) a typical random configuration of the phases {φ n }, (ii) the case of the asymmetric advection field (3) with no disorder φ n = 0 ∀n, iii) the case of the symmetric advection field (4) with no disorder φ n = 0 ∀n.
The four cases analyzed show the presence of three distinct scaling regimes for t D . In the case of average over the disorder, the moments [x(t) − x(0)] m w φ display a typical scaling form of the type ∝ t αm in the three regimes, only for even values of m. On the other hand, they are equal to 0 for odd m. The scaling of the moments [x(t) − x(0)] m w 0 when a generic random configuration of the phases {φ n } has been set up, traces that of [x(t) − x(0)] m ] w φ only when m is even. The odd moments show random sample-to-sample fluctuations which, obviously, are not consistent with any scaling law. In this sense, we can state that the average over the initial conditions for a generic quenched disorder is typical, and follows a universal scaling. The cases of zero disorder exhibit scaling trends differ from those of a generic random configuration for times t L 2 M 2 D . Specifically, the symmetric advection field (4) has moments whose scaling form is strongly anomalous for any m. In the antisymmetric case (3), the even moments behave exactly like the symmetric counterparts. The odd moments, on the other side, are random functions of time for which a clear scaling form cannot be drawn.
Finally, in this article, we offered the most detailed study of the scaling form of the position moments in the random shear model. Our findings include and extend to all the previous results present in the literature, providing a comprehensive and accurate theoretical framework on this subject. Moreover, it constitutes a benchmark for future analysis including the analytical calculations of the propagator P(x, t), the diffusion equation that it satisfies and the question of the ergodicity breaking in either one of the two averaging procedure.

Conflicts of Interest:
The authors declare no conflict of interest.