The Eulerian space-time correlation of strong Magnetohydrodynamic Turbulence

The Eulerian space-time correlation of strong Magnetohydrodynamic (MHD) turbulence in strongly magnetized plasmas is investigated by means of direct numerical simulations of Reduced MHD turbulence and phenomenological modeling. Two new important results follow from the simulations: 1) counter-propagating Alfv\'enic fluctuations at a each scale decorrelate in time at the same rate in both balanced and imbalanced turbulence; and 2) the scaling with wavenumber of the decorrelation rate is consistent with pure hydrodynamic sweeping of small-scale structures by the fluctuating velocity of the energy-containing scales. An explanation of the simulation results is proposed in the context of a recent phenomenological MHD model introduced by Bourouaine and Perez 2019 (BP19) when restricted to the strong turbulence regime. The model predicts that the two-time power spectrum exhibits an universal, self-similar behavior that is solely determined by the probability distribution function of random velocities in the energy-containing range. Understanding the scale-dependent temporal evolution of the space-time turbulence correlation as well as its associated universal properties is essential in the analysis and interpretation of spacecraft observations, such as the recently launched Parker Solar Probe (PSP).


I. INTRODUCTION
Since the first observation that non-compressive Alfvénlike fluctuations of velocity and magnetic field dominate the solar wind by Belcher and Davis [1], incompressible Magnetohydrodynamics (MHD) [2] has been often invoked to describe the observed Kolmogorov-like the power spectrum of low-frequency fluctuations of the solar wind plasma, for an extensive review see Bruno and Carbone [3], Chen [4], Verscharen et al. [5]. The majority of advances in MHD turbulence in the last few decades have been largely concerned with its spatial statistical properties, such as the three dimensional structure of the power spectrum and higher order structure functions [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Most of these properties can be derived from two-point one-time correlations, which quantify the covariance between simultaneous values of a turbulent quantity at two different points.
However, more often than not turbulence experiments and solar wind observations can only provide single-probe measurements along the plasma at different times and locations, in which case a methodology to relate the time signals measured in the probe-frame to the spatial properties in the plasmaframe is required to test theoretical predictions. For instance, in solar wind observations the so called Taylor Hypothesis (TH) [20] (or frozen-in-flow approximation) is commonly used. This approximation essentially establishes that when the mean flow speed U (as seen in the probe frame) is much larger than any other characteristic speed, such as the flow's turbulent amplitude u 0 and wave-propagation speed, the time signal of a turbulent quantity measured in the probe's frame is due to the advection of a frozen spatial structure passing by the instrument at the local flow speed U. In contrast, when these conditions are violated the temporal variation observed singleprobe measurements arise instead from a dynamic structure passing by the probe, i.e., the time variation is a combination of advection and evolution of the passing structures. The recently launched Parker Solar Probe (PSP) [21] has spurred a renewed interest in understanding the space-time structure of solar wind turbulence [22][23][24][25][26], precisely because it will explore the near-Sun region where the conditions for the validity of the TH might not be satisfied. In this case, an understanding of the structure of two-time two-point correlations of turbulent quantities and any possible universal properties are essential to successfully relate the turbulent time signals (measured by the probe) to the spatial structure of the turbulence. Analyses of spacecraft data to date have provided increasing evidence that many turbulent properties of low frequency fluctuations in the solar wind are consistent with various predictions of current MHD turbulence models [4], however the subject remains open and under active investigation and debate. In this paper we address the problem of the physics of temporal decorrelation of Alfvénic fluctuations, which may be relevant to the analysis of solar wind fluctuations whenever they can be described by incompressible MHD [22][23][24][25][26][27][28], such as in the first two perihelia around 36 solar radii where the solar wind was found to be highly Alfvénic [29][30][31][32][33][34][35]. The physics of the temporal decorrelation in solar wind observations when MHD is not applicable is outside the scope of this work and deserves further investigation.
A number of works on the structure of the Eulerian spacetime correlation in MHD turbulence have been carried out in the past decade. Servidio et al. [27] investigated the scaledependent temporal correlation in isotropic MHD using numerical simulations. The authors reported that the Eulerian decorrelation time is consistent with a sweeping-like scaling τ c ∼ k −1 , which they attribute to a combination of convective sweeping by the large-scale flow and magnetic sweeping from large-scale magnetic fluctuations. Lugones et al. [28] studied the same problem for the anisotropic case of MHD turbulence with a guide field, using simulations with small, moderate and large guide field to investigate the role of the magnetic field in the decorrelation time. Their findings are consistent with Servidio et al. [27] for small fields, while the decorrelation becomes dominated by Alfvén-wave propagation for large guide field. One shortcoming of the works of Servidio et al. [27] and Lugones et al. [28] is that they focused on the correlation function of the fluctuating magnetic field and not on the Elsasser fields. Narita [36] investigated the temporal decorrelations of the Elsasser fields in MHD, by extending a Hydrodynamic (HD) sweeping model of the Eulerian correlation of Wilczek and Narita [37] with a mean flow. His model suggests that the temporal decorrelations of the Elsasser fields z ± at small-scales arises from the random sweeping by large-scale Elsasser fluctuations propagating in the opposite direction z ∓ , resulting in different decorrelation rates for imbalanced turbulence. However, Bourouaine and Perez [23] measured the Eulerian correlation of Elsasser fields in highly imbalanced, reflection-driven MHD turbulence simulations with high space-time resolution and found that the decorrelation of both fields is consistent with sweeping by large-scale fluctuations at a common speed that is comparable to the root mean squared (r.m.s.) value of the fluctuating velocity, suggesting that the sweeping is hydrodynamic in nature. Based on the evidence from the numerical simulations, Bourouaine and Perez [24, hereafter BP19] recently introduced a new sweeping model of MHD turbulence that relies on the local mean field direction. The findings from the model are consistent with a common sweeping characteristic timescale for both Elsasser fields. In this work we show that when the BP19 model is applied in strong MHD turbulence, the Eulerian space-time correlation is entirely dominated by HD sweeping, which we confirm in numerical simulations.
This paper is organized as follows. In section II we discuss the theoretical framework for the simulations and phenomenological models presented in this work, including a brief description of Kraichnan's idealized convection model in HD and its generalization to strong MHD turbulence. In section III we discuss the numerical simulation setup and simulation parameters and a brief description of the methodology that will be used to validate the MHD sweeping model in simulations. In section IV we present and discuss the simulation results and in section V we conclude.

II. THEORETICAL FRAMEWORK
We assume that velocity v(x,t) and magnetic field B(x,t) fluctuations are described by the equations of ideal incompressible MHD, which in terms of the Elsasser variables z ± ≡ v ± b take the form where b = B/(4πρ) is the fluctuating Alfvén velocity, V A = B 0 /(4πρ) is the background Alfvén velocity, ρ is the constant background plasma density and p is the combined thermal and magnetic pressure. Random forcing f ± and viscous dissipation terms have been included to investigate the case of steadily-driven turbulence. For a strong background magnetic field (B 0 = |B 0 |ê z , say, with |B 0 | ≫ |b|) the universal properties of MHD turbulence can be accurately described by neglecting the field-parallel component, z ± , of the fluctuating fields (the pseudo-Alfvén fluctuations) that play a sub-dominant role in the turbulence dynamics (see [13] and references therein). It can be further demonstrated that setting z = 0 in equation (1) leads to a set of equations that is equivalent the simpler Reduced MHD (RMHD) model [38,39]. It is worth noting that the RMHD model is commonly invoked to describe the dominant nonlinear interactions and resulting turbulence of noncompressive Alfvén-like fluctuations, which comprise most of the energy in the solar wind. It has also been shown, from gyrokinetics [40] and from comparisons with MHD simulations [14,41], that RMHD rigorously describes the essential non-linear interactions responsible for the turbulence cascade of non-compressive Alfvénic fluctuations.
For homomogeneous and stationary Elsasser fluctuations z ± (x,t) the two-time two-point Eulerian correlation is defined as where · · · represents an ensemble average over many turbulence realizations. These correlations measure the degree to which each Elsasser field at any position x and time t is correlated with itself at another location with relative position r after a time τ has elapsed. In a turbulent system correlations arise due to the presence of coherent structures of many characteristic length-scales, which are undergoing random advection and nonlinear straining by other structures in the flow according to the dynamics determined by equations (1). Although turbulence correlations can, at least formally, be related to the governing equations of the fluctuating variables, turbulence theories have been unable to produce exact analytical solutions even in the simplest case of incompressible HD turbulence, because of the well known closure problem [42]. The correlation in equation (2) can be expressed in terms of its spatial Fourier transform where h ± (k, τ) are the so-called two-time power spectra and z ± (k,t) is the spatial Fourier transform of the field z ± (x,t). The scaling properties and three-dimensional structure of the spatial power spectra h ± 0 (k) = h ± (k, 0) (for τ = 0) have been the subject of extensive research in theory, numerical simulations and solar wind observations [4]. In this work, we make very few assumptions about the structure of the spatial power spectrum and focus our investigation on the structure of the scale-by-scale τ dependency, which accounts for the scale-dependent temporal decorrelation of the turbulence. As τ increases, the Fourier amplitudes at wavevector k decorrelate and one can thus define the scale-dependent time correlations Γ ± (k, τ) as By definition Γ ± (k, 0) = 1 for all k, which means that at zero time lag τ the fluctuations are perfectly correlated. The advantage of equation (5) is that it allows for the separation of the spatial part of the correlation function from the scaledependent temporal part, which is the subject of this work. The two-time power spectra h ± (k, τ) are simply a different representation of the correlation functions C ± (r, τ) and thus contain the same information.

A. Kraichnan's idealized convection model in Hydrodynamics
In HD turbulence, temporal decorrelation at a given point can arise from two main effects: 1) random sweeping of the small-scale eddies by large ones and 2) eddy straining (or shear) associated with nonlinear inertial forces. Scaling arguments can be used to argue that the Eulerian correlation in Hydrodynamic (HD) is dominated by the first of these two effects, also known as the Kraichnan's Sweeping Hypothesis (KSH) [43]. The sweeping decorrelation mechanism is a non-local-in-scale process, in the sense that it involves eddies of disparate scales, and its characteristic timescale is τ s ∼ 1/(ku 0 ) for a fluctuation of scale λ ∼ 1/k swept by a large-scale fluctuation with velocity u 0 . The second timescale associated with nonlinear straining scales as τ NL ∼ k −2/3 , a slower decrease with k than the sweeping timescale τ s . These scaling arguments suggest that the Kraichnan's hypothesis is expected to hold better for sufficiently large values of k for which the ratio τ s /τ NL ∼ k −1/3 is small.
Kraichnan introduced an idealized convection model of incompressible HD to describe the random sweeping of smallscale fluctuations by large ones. In this model the fluid velocity consists of two parts v = v ′ + u, with the following assumptions: 1) v ′ , describing the large-scale eddies, is constant in space and time but is a zero-mean random variable with an isotropic Gaussian distribution, 2) the field u(x,t), describing the small-scale eddies, is much smaller in magnitude than v ′ , and 3) v ′ and u(x, 0) are statistically independent. From the first two assumptions, the Navier-Stokes equation becomes where we ignored viscous dissipation (considering fluctuations in the inertial range) and dropped the pressure term whose only role is to ensure fluctuations remain incompressible. As opposed to the Navier-Stokes (NS) equation, the idealized model given by equation (6) is a stochastic linear equation in u, which does not have the statistical closure problem of the nonlinear NS equation. The essence of this approximated idealized model is that the dominant variation of u simply arises from advection of frozen structures by a constant but random velocity at each point. For instance, if we assume for the moment that v ′ = V is not a random variable, equation (6) forms the basis for the TH approximation. In this sense, the random sweeping model can be interpreted as the application of TH to a statistical ensemble of systems, each one with a different large scale flow velocity drawn from a random distribution corresponding to the large-scale eddies. Straightforward solutions to Equation (6) can be found to obtain the scale-dependent time correlation where P(v ′ ) is the probability density for the random variable v ′ . Kraichnan's model assumed P(v ′ ) to be an isotropic Gaussian distribution in which case equation (7) becomes where γ k ≡ kv 0 / √ 2 is the decorrelation rate and v 0 is the r.m.s. value of the velocity v ′ along any given direction. The decorrelation rate is defined at each k as γ k = 1/τ k , where τ k is the time lag for which the correlation Γ(k, τ) drops to 1/e ≃ 0.37. This idealized model provides a phenomenological description of the temporal decorrelation when the timescale τ c ∼ 1/(kv 0 ) is much faster than the Kolmogorov estimate of the nonlinear cascade time τ NL ∼ k −2/3 . Wilczek and Narita [37] revisited the HD case with a constant mean flow U, which simply adds a phase factor to the correlation Γ(k, τ).
It is important to note that Kraichnan's assumption of Gaussianity for the random variable v ′ is not necessary and the validity of his model can be extended to other distributions P(v ′ ) by noticing that the average in equation (7) is nothing but the characteristic function ϕ v ′ (ξ) of the probability density P(v ′ ), hence where This result shows that the scale-dependent time correlation can be obtained from characteristic function of the probability density of the large-scale eddies [44], by setting the velocitywavenumber ξ equal to kτ, and is therefore self-similar. In the next subsection we extend this idealized model for the case of strong MHD turbulence following the phenomenology introduced by BP19. B. Sweeping model for strong MHD turbulence.
The Kraichnan's picture acquires greater complexity in MHD turbulence for a number of reasons. First, in the Elsasser formulation MHD contains two fluctuating fields z ± that are being advected in opposite directions along the background magnetic field and undergo mutual straining only when counter-propagating fields encounter each other or "collide", resulting in various limiting regimes. For instance, when the energies of the fluctuating fields z + and z − are comparable the turbulence is called balanced, otherwise it is called imbalanced. For both the balanced and imbalanced cases, the turbulence can be weak [6] or strong [45]. The weak turbulence regime occurs when the time it takes two eddies to cross one another is much shorter than the nonlinear interaction time, thereby requiring a large number of successive collisions before eddies can cascade their energy to smaller scales [6]. In the strong regime, the crossing and nonlinear times are comparable and the cascade occurs in a single collision. Although it is still a matter of debate, a number of models assume that for imbalanced turbulence z ± may have Run Regime Normalized cross-helicity (σ c ) u 0 b 0 z + 0 z different nonlinear straining times τ ± NL ∼ λ /z ∓ λ , due to the amplitude difference [8,46,47].
BP19 introduced a new model for the scale-dependent time correlations, by writing the large-scale Elsasser variables z ′ ± = v ′ ± b ′ in terms of their corresponding fluctuating velocity and magnetic fields in (1) to obtain where v ′ , b ′ and δ z ± are taken to be perpendicular to the magnetic field, consistent with the RMHD approximation. Here δ z ± represent the small-scale Elsasser fluctuations, V A ′ ≡ V A + b ′ is the modified Alfvén velocity resulting from the superposition of the mean background magnetic field and the fluctuating component b ′ from the large-scale eddies, and ∇ ⊥ , ∇ are the field-perpendicular and field-parallel gradient operator defined with respect to the local magnetic field, respectively. Hereafter, primes are used to represent random variables with known statistics. Under the assumption that the characteristic timescales of the RHS terms are much smaller than those in the LHS, for a strongly magnetized plasma (|b ′ | ≪ V A ) and for Gaussian-distributed outer-scale velocities v ′ one obtains [24] Γ ± (k, τ) = e ∓ik V A τ e −(γ k τ) 2 (13) where V A is the Alfvén speed, k is the component of k in the direction of the local magnetic field and γ k = k ⊥ v 0 / √ 2 is the decorrelation rate. This result is very similar to the model obtained by Narita [36], with the important difference that both Elsasser fields decorrelate at a common rate, determined by pure HD sweeping. Noting that v 0 represents the r.m.s. of the fluctuating velocity in any direction and v ′ lies in the fieldperpendicular plane, the velocity r.m.s. is u 0 = √ 2v 0 in which case γ k = k ⊥ u 0 /2.
Scaling arguments can also be used to obtain a model for the Γ ± (k, τ) functions in the strong turbulence regime. Let λ and l be the field-perpendicular and field-parallel lengthscales of an eddy with respect to the local magnetic field and of amplitude v λ . If the turbulence is driven isotropically, λ ∼ l, the turbulence is necessarily weak when v λ ≪ V A . Because weak turbulence cascades energy to smaller perpendicular scales without affecting the parallel structure l [6], eddies will progressively become elongated along the field until the nonlinear time τ λ ∼ λ /v λ becomes comparable to the linear timescale τ A ∼ l/V A . Therefore, the turbulence will unavoidably become strong when the critical balance condition [45] τ λ ∼ τ A ⇒ λ /v λ ∼ l/V A is satisfied, which means the timescale of Alfvén wave propagation becomes comparable to that of the nonlinear terms in the RHS of equation (12). In this case, the Alfvénic propagation can be neglected and the correlation function takes the form and in the case of a Gaussian distribution P(v ′ ) Equations (14) and (15) show that the decorrelation is solely determined by HD sweeping and is the same for both Elsasser fields.
It is worth mentioning here that the assumption of a Gaussian distribution of large-scale velocities is made here for concreteness, to make calculations simpler and because the simulations used in next section to validate the model are driven at the outer-scale in a Gaussian fashion. However, for the application of this model to solar wind observations [see for instance 26], the actual distribution of velocities at the outerscale can be used. Previous solar wind observations have shown evidence that large-scale velocities and magnetic field fluctuations are Gaussian [3], although some observations suggest these fluctuations can show strongly non-Gaussian, skewed tails [48].

III. NUMERICAL SIMULATIONS
RMHD equations, obtained by setting z = 0 in Eq. (1), are solved using a fully dealiased 3D pseudo-spectral code in a rectangular domain with aspect ratio 1 : M, where M is the ratio between parallel and perpendicular box sizes, defined with respect to the background magnetic field B 0 = B 0êz . The normalization chosen in the simulations is such that the r.m.s. values of fluctuating plasma and Alfvén velocity are of order u rms ∼ 1 (in code's units), and the magnitude B 0 so that V A /u rms is of order M. The box size in the xy plane (perpendicular to the guide magnetic field) is chosen as L = 2π and time is normalized to the large scale eddy turnover time τ 0 = L/2πu rms . The turbulence is driven by random forcing in the field-perpendicular and field-parallel wave-numbers 0 < k ⊥ < 4 and 0 < k < 2, respectively, which due to the aspect ratio of the simulation box allows one to drive strong RMHD turbulence by controlling the degree to which outer-scale eddies satisfy the critical balance condition k V A ∼ k ⊥ u rms [49], where u rms represents the r.m.s. value of the turbulent velocity. The Reynolds number is defined as Re = u rms (L/2π)/ν.
Approximately every eddy turnover time in the steady state the code outputs snapshots of the Elsasser fields z ± , as well as the entire time history of each Elsasser field on eight selected xy planes. In the simulations, correlations between v and b are introduced through the random forcing to investigate the role of cross-helicity H c = E + − E − , which measures the energy difference between counter-propagating Elsasser fluctuations. Cross-helicity is conveniently quantified in the simulations through the normalized cross-helicity σ c = H c /E defined as the amount of cross-helicity normalized to the total energy E = E + + E − . The normalized cross-helicity takes values in the range −1 ≤ σ c ≤ 1, with zero corresponding to balanced turbulence, and imbalanced turbulence otherwise.
Three simulations of steadily-driven RMHD turbulence, listed in table I, are used to investigate the scaling properties of the time correlations Γ ± (k, τ) and to compare with phenomenological models . These simulations have been extensively used to investigate the structure and scaling of the spatial spectrum h 0 (k) of balanced and imbalanced MHD turbulence in previous works [11][12][13]. Because the parameters of the simulations in table I are the same as simulations RB1, RB2 and RI1 in Perez et al. [13], we adopt the same labeling convention.
The random forcing drives the outer-scale velocities toward an isotropic two-dimensional Gaussian distribution of the form given in equation (8), and whose characteristic function is where u 0 is the r.m.s. value of velocities in the energycontaining range and ξ is a velocity-wavevector on which the characteristic function depends. The Gaussian nature of the outer-scale flow is verified by measuring the angle-averaged characteristic function ϕ v ′ (ξ ) = ϕ v ′ (ξ) φ for the most energetic fluctuations, taken as wavenumbers k ⊥ 7. Figure 1 shows that when the characteristic functions ϕ v ′ for all three simulations are plotted versusξ ≡ ξ u 0 , they all overlap almost exactly with the Gaussian function represented by the circles in the figure.
The two-dimensional, two-time power spectrum of each Elsasser field is calculated through the average in terms of the Fourier transforms z ± α ≡ z ± (k ⊥ , ζ α ,t) at each transverse plane ζ α , with α = 1, · · · 8. The ensemble average · · · φ in this equation also includes an average over the polar angle in the k ⊥ plane, due to the isotropy of the two-dimensional power in the field-perpendicular plane. The scale-dependent time correlations can be related to the  two-dimensional two-time power spectra by integrating equation (5) over k Assuming that the scale-dependent time correlations Γ ± (k ⊥ , k , τ) weakly depend on k as predicted by the model it follows that where is the two-dimensional power spectrum. The fieldperpendicular energy spectra, defined as measured for simulations RB1, RB2 and RI1 are shown in Figure 2. As previously reported, the simulations are consistent with the scale-dependent phenomenology of strong MHD turbulence [9] for balanced and imbalanced turbulence [12].   In the next section we present results from the scaledependent time correlations measured in the simulations from equation (21), and perform two important tests to compare with the theoretical formula given in equations (14) and (15) . The first test is to show that decorrelation rates γ ± k scale linearly with k ⊥ , consistent with sweeping characterized by the r.m.s. speed of the outer-scale flow. In the second test the scale-dependent time correlations Γ ± (k ⊥ , τ), as defined in equation (21), are calculated in numerical simulations for wavenumbers k ⊥ in the inertial range. The resulting correlations are then compared with the characteristic function ϕ v ′ (k ⊥ τ), computed from equation (11) using the random distribution of outer-scale velocities P(v ′ ).

IV. SIMULATION RESULTS.
The modeled scale-dependent time correlations for strong MHD turbulence given in equation 15 have the following important features: 1) they are solely due to the random sweeping by the large-scale flow, 2) the decorrelation rates are the same for both Elsasser fluctuations z ± , whether the turbulence is balanced or imbalanced, and scale linearly with k ⊥ as γ ± = k ⊥ u 0 /2 in the inertial range, and 3) they exhibit universal, self-similar behavior as they can all be written in terms  Figure 4. Sweeping velocities c ± = 2γ ± k /k ⊥ estimated from measured decorrelation rates γ ± k vs k ⊥ , obtained from Gaussian fits shown in Figure 3. The decorrelation rates are similar for both Elsasser fields at all wavenumbers k ⊥ , for both balanced and imbalanced turbulence. Flat regions of c ± correspond to wavenumbers for which the decorrelation is consistent with sweeping, and the char- The scale-dependent time correlations Γ ± (k ⊥ , τ) are calculated for each simulation according to equation (21) and shown in Figure 3 for selected values of the fieldperpendicular wavenumber between k ⊥ = 16 and k ⊥ = 64. Dotted lines correspond to Gaussian fits of each correlation of the form given in Equation (15) where the only free parameter is the corresponding decorrelation rate γ ± k . As noted by Bourouaine and Perez [23] in simulations of reflectiondriven Alfvén turbulence, the scale-dependent time correlations closely follow Gaussian behavior in τ, consistent with the statistics associated with the energy containing scales. Figure 4 shows c ± ≡ 2γ ± k /k ⊥ vs k ⊥ from the measured decorrelation rates for each Elsasser field z ± , which is consistent with the prediction of BP19 in the strong turbulence regime in a number of ways. First, the decorrelation rates for both fields remain approximately equal at all wavenumbers for both balanced and imbalanced simulations. Second, c ± approximately starts to plateau at wavenumbers above k ⊥ 4 and remains approximately constant up to wavenum-  bers around k ⊥ ≃ 20 for simulations RB1 & RI1 and up to k ⊥ ≃ 30 for RB2, which indicates the decorrelation rates exhibit linear behavior consistent with sweeping in the inertialrange scales for each simulation. Moreover, the theoretical model predicts that c ± is a common speed for both Elsasser variables z + and z − and equal to the characteristic sweeping speed u 0 . Figure 4 shows four horizontal lines corresponding to the r.m.s. values of v, b, z ± of the most energetic eddies, and c ± is clearly consistent with the r.m.s. of the outer-scale velocities only, as indicated by the solid horizontal line.
The second test of the theoretical model is found in Figure 5, which shows the computed scale-dependent time correlation functions for selected values of k ⊥ in runs RB2 & RI1 of balanced and imbalanced RMHD turbulence. In all these plots it is observed that when the corresponding temporal correlation function is plotted vs the normalized velocitywavenumberξ = k ⊥ τu 0 , the scale dependent correlations for all wavenumbers are essentially indistinguishable from the characteristic function arising from the random distribution of velocities of the most energetic scales, as predicted by the model in equation (14).

V. CONCLUSION.
In this work we investigated the temporal decorrelation of strong MHD turbulence through phenomenological modeling and numerical simulations of RMHD turbulence. Scaledependent time correlations of Elsasser fluctuations were modeled by restricting the recent BP19 phenomenology to the strong MHD turbulence regime. In the BP19 phenomenol-ogy, which extends Kraichnan's idealized convection model of Hydrodynamics to MHD, the Eulerian decorrelation is the result of HD sweeping and Alfvénic propagation along the local magnetic field. In this work we have shown that, for the particular case of strong turbulence regime, the decorrelation is solely dominated by HD sweeping. The resulting scaledependent time correlations exhibit an universal, self-similar behavior that is entirely determined by the statistics of velocities at the largest, energy-containing eddies, and it clearly shows that the decorrelation rates for both Elssaser variables are the same regardless of whether the turbulence is balanced or imbalanced. All these features were tested using numerical simulations of strong RMHD turbulence. The numerical results are in very good agreement with the theoretical predictions.
An earlier model of the Eulerian decorrelation in MHD was also proposed by Narita [36] in which the decorrelation rates are predicted to scale as γ ± k = k ⊥ z ∓ rms / √ 2. This model is physically appealing at a first glance, as it is natural to assume that the decorrelation rate of z + is determined by the r.m.s. value of z ′ − , and viceversa. However, as [36] points out, the decorrelations rates will be different for imbalanced turbulence where one of the Elsasser component has a larger amplitude than the other. One shortcoming of Narita [36] model is that it does not take into account the fundamentally different effects that the large-scale flow and magnetic field fluctuations have on the small-scale ones. For instance, large-scale fluctuations in velocity will sweep small-scale eddies equally for both Elsasser fields, while large-scale fluctuations in the magnetic field simply modify the background magnetic field along which smallscale eddies propagate. It is therefore important that a sweeping model captures these different characteristics associated with the nature of large-scale velocity and magnetic field fluctuations.
The outcomes from this study will in fact be very beneficial for the analysis of the spacecraft signals beyond the validity of Taylor's Hypothesis whenever solar wind observations are compared with predictions from phenomenological MHD turbulence models [see for instance 26]. Finally we conjecture that the results we present in this work may be suitably extended to kinetic (non-MHD) regimes found in the solar wind, but it requires further investigation.