4.1 Turbulent-decay phenomenology

The basic idea of the phenomenological model is to treat the dominant $\tilde {z}^+$ fluctuations as a standard decaying turbulence problem, while $\tilde {z}^{-}$ is effectively strongly forced by reflection and damped by turbulence. In more detail, because $\tilde {E}^+ \gg \tilde {E}^{r}$ when $\tilde {E}^+\gg \tilde {E}^{-}$ (as assumed), reflection is negligible for the $\tilde {z}^+$ field, and consequently, for the forcing/damping of the wave-action energy (see (2.8)). This implies that $\tilde {E}^+$ is approximately ideally conserved during this phase and its turbulent decay occurs only due to nonlinear coupling with $\tilde {z}^{-}$. Throughout this phase, the $\tilde {z}^{-}$ fluctuations, which are forced by reflections, remain very low amplitude; therefore a-priori, one might expect $\tilde {z}^+$ fluctuations to be in the weak regime. However, we assume (Velli et al. Reference Velli, Grappin and Mangeney1989; Lithwick, Goldreich & Sridhar Reference Lithwick, Goldreich and Sridhar2007; Chandran & Perez Reference Chandran and Perez2019; Schekochihin Reference Schekochihin2022), providing detailed numerical justification below (§ 4.3), that the $\tilde {z}^{-}$ fluctuations remain ‘anomalously coherent’ with the $\tilde {z}^+$, because their forcing via reflection is highly coherent (${\propto }-\tilde {\boldsymbol {z}}^+$), thus ‘dragging’ $\tilde {z}^{-}$ along with the $\tilde {z}^+$ in time. The consequence is twofold: first, by moving into the frame that propagates outwards with $\tilde {\boldsymbol {z}}^+$, it allows one to ignore the Alfvénic propagation terms for both $\tilde {\boldsymbol {z}}^+$ and $\tilde {\boldsymbol {z}}^{-}$; second, it allows the estimation of turbulent cascade times using the standard nonlinear turnover times (unlike for weak turbulence). Therefore, the turbulent-decay time $\tau ^{\pm }$ of $\tilde {z}^{\pm }$ is

(4.1)\begin{equation} \tau_{{\mp}}^{{-}1}\sim a^{{-}3/2}\frac{\tilde{z}^{{\pm}}}{\tilde{\lambda}^{{\pm}}} = \frac{z^{{\pm}}}{\lambda^{{\pm}}}, \end{equation}

where $\tilde {\lambda }^{\pm }$ are the characteristic perpendicular scales of $\tilde {z}^{\pm }$ in the co-moving frame that govern the decay/growth of $\tilde {z}^{\mp }$ and $\tilde {z}^{\pm }$ represents the r.m.s.amplitude of $\tilde {\boldsymbol {z}}^{\pm }$. Variables without the tilde are in the physical frame with physical units, showing how the $a^{-3/2}$ factor arises from the use of wave-action variables.

Based on these assumptions, we compute the evolution of $\tilde {z}^{-}$ via the balance of reflection and nonlinear decay, ignoring the Alfvénic and time-evolution terms (the latter is small, as justified below). The evolution of $\tilde {z}^+$ results from its nonlinear turbulent decay via the $\tilde {z}^{-}$ that it has sourced. The scheme then yields the following phenomenological evolution equations for $\tilde {E}^{\pm }$ (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002):

(4.2a)\begin{gather} \dot{a}\dfrac{\partial {\tilde{E}^+}}{\partial a} \sim{-}\dfrac{1}{a^{3/2}}\dfrac{\tilde{z}^{-}}{\tilde{\lambda}_{-}} \tilde{E}^+, \end{gather}

(4.2b)\begin{gather}\dfrac{1}{a^{3/2}}\dfrac{\tilde{z}^+}{\tilde{\lambda}_+}\tilde{E}^{-}\sim \dfrac{\dot{a}}{a}|\tilde{E}^{r}|\sim \frac{\dot{a}}{a}|\sigma_{\theta}| \tilde{z}^+ \tilde{z}^{-}. \end{gather}

Writing (4.2b) for $\tilde {z}^{-}$ instead gives

(4.3)\begin{equation} \tilde{z}^{-} \sim \dot{a}a^{1/2} \tilde{\lambda}_+|\sigma_{\theta}|, \end{equation}

whereby we see the interesting feature that the amplitude of $\tilde {z}^{-}$ is independent of that of $\tilde {z}^+$ (other than indirectly through $\tilde {\lambda }_+$ and $\sigma _{\theta }$). This occurs because $\tilde {z}^+$ acts to both drive and dissipate the $\tilde {z}^{-}$ energy. This independence from the $\tilde {z}^+$ spectrum also suggests that, with various caveats discussed below (§ 4.2), it could be reasonable to reinterpret the balance of reflection and nonlinear damping as applying at each scale separately, thus replacing the $\tilde {\lambda }_+$ in (4.3) with $\tilde {k}_{\perp }^{-1}$ and making $\tilde {z}^{-}$ the r.m.s. amplitude of the $\tilde {z}^{-}$ increment across a distance $\tilde {k}_\perp ^{-1}$ in the perpendicular plane. This gives $\tilde {z}^{-}(k_{\perp })\propto k_{\perp }^{-1}$ or a ${\propto } k_{\perp }^{-3}$ energy spectrum for $\tilde {z}^{-}$. We can insert (4.3) into (4.2a) to obtain the total energy ($\tilde {E}\approx \tilde {E}^+$) decay,

(4.4)\begin{equation} \frac{\partial \ln \tilde{E}^+}{\partial a} \sim{-}\frac{1}{a} \frac{\tilde{\lambda}_+}{\tilde{\lambda}_{-}}\sigma_{\theta}.\end{equation}

Several other points are worth noting. First, the anomalous coherence will break down once the effect of $\tilde {z}^+$ on $\tilde {z}^{-}$ enters the weak regime (in which case $\tilde {z}^{-}$ can propagate away from its $\tilde {z}^+$ source). The phenomenology thus requires

(4.5)\begin{equation} \chi_{A}\doteq \frac{(\tau_{-})^{{-}1}}{ v_{A}/\ell_{\|}} \sim \frac{\tilde{z}^+{/}\tilde{\lambda}_+}{a^{1/2}v_{A0}/\ell_{\|}}\gtrsim1, \end{equation}

where $\ell _{\|}$ is the parallel correlation length ($\chi _{A}>1$ may be unphysical for other reasons, but the phenomenology itself is fundamentally two dimensional, ignoring $\ell _{\|}$). Second, we verify that the neglect of $\partial _{t} \tilde {z}^{-}$ is consistent, so long as anomalous coherence allows us to ignore the Alfvénic propagation of $\tilde {z}^{-}$ in the frame of $\tilde {z}^+$, by noting that $\partial _{t} \tilde {z}^{-} \sim (\dot {a}/a)\tilde {z}^{-}$ is a factor ${\sim }\tilde {z}^{-}/\tilde {z}^+$ smaller than the reflection term in (4.2b). Third, there exists an additional constraint implicit in (4.2), which comes from noting that (4.3) is equivalent to

(4.6)\begin{equation} \tilde{z}^{-} \sim \frac{\tilde{z}^+}{\chi_{\rm exp}},\end{equation}

where

(4.7)\begin{equation} \chi_{\rm exp}\doteq \frac{(\tau_{-})^{{-}1}}{\dot{a}/a} \sim \frac{\tilde{z}^+{/}\tilde{\lambda}_+}{a^{1/2}\dot{a}}\end{equation}

is the ratio of the nonlinear damping to reflection rates. Thus, the phenomenology can only be valid for sufficiently large-amplitude $\tilde {z}^+$ with $\chi _{\rm exp}\gg 1$, irrespective of the fluctuation's parallel scale, and we expect the transition to the balanced regime to occur when $\tilde {E}^+$ decays sufficiently so that $\chi _{\rm exp}\sim 1$. Here $\chi _{\rm exp}$ will feature prominently below as the key parameter that controls the transition from the imbalanced to balanced phase.

Previous treatments (Verdini & Velli Reference Verdini and Velli2007; Chandran & Hollweg Reference Chandran and Hollweg2009; Verdini et al. Reference Verdini, Velli, Matthaeus, Oughton and Dmitruk2010) have taken $\tilde {\lambda }_+$ and $\tilde {\lambda }_{-}$ in (4.2) to be the same and constant in time in the co-moving frame. But, the argument about the $\tilde {z}^{-}$ balance and spectrum in the previous paragraph, as well as decaying turbulence theory in general (Kolmogorov Reference Kolmogorov1941), suggest that there is no reason to expect this to be the case. Indeed, if the $\tilde {z}^{-}$ spectrum was ${\propto }k_{\perp }^{-3}$ as suggested above, then – irrespective of the dominant scales of $\tilde {z}^+$ – the correlation scale of $\tilde {z}^{-}$ would become the largest scale at which the arguments leading to (4.2) break down (e.g. where $\chi _{\rm exp}<1$, or where the turbulence becomes weak). In addition, we will show below that the co-moving scales of $\tilde {z}^+$ evolve in time as a result of another nonlinear conservation law (that for the wave-action anastrophy). Herein lies the problem that complicates the comparison of the phenomenology to the numerical experiments: it is not clear what constrains the $\tilde {\lambda }_+$ and $\tilde {\lambda }_{-}$ scales in (4.2), but their time evolution is crucial for determining many key aspects of the turbulent evolution. In addition, it is not clear how the evolution of $\tilde {\lambda }_+$, which is the characteristic scale of $\tilde {z}^+$ that controls the nonlinear evolution of $\tilde {z}^{-}$, relates to that of the correlation scale $\tilde {L}_+$ of $\tilde {z}^+$. This allows us to consider the evolution of $\tilde {L}_+$ separately from the decay phenomenology, unlike in standard decaying turbulence theory (§ 4.4), but the cause of this apparent discrepancy between $\tilde {\lambda }_+$ and $\tilde {L}_+$ remains a poorly understood aspect of the phenomenology.

4.2 Turbulent spectra

The energy spectra $\tilde {\mathcal {E}}^{\pm }(k_{\perp })$ for the $\chi _{A0}=96$ ($\chi _{\rm exp0}=960$) simulation over this imbalanced phase are shown in figure 3. Panel (a) shows the time evolution of $\tilde {\mathcal {E}}^+$ and $\tilde {\mathcal {E}}^{-}$, respectively, demonstrating their very different evolution. Panel (b) shows the simulation at $a=5.2$ when it has been refined to a resolution $n_{\perp }^{2}\times n_{z}=8192^{2}\times 256$ in order to attempt to capture the transition to standard imbalanced turbulence at small scales. Footnote ⁷ The obvious feature of $\tilde {\mathcal {E}}^+(k_{\perp })$ is its rapid migration towards large scales, which will be discussed in detail below in § 4.4. As this occurs, $\tilde {\mathcal {E}}^+$ develops a wide $\tilde {\mathcal {E}}^+\propto k_{\perp }^{-1}$ range, which eventually transitions into a steeper slope at small scales (see panel (c) at $a\approx 5.2$). While the simulation does not have sufficient resolution to easily diagnose the slope of this smaller-scale turbulence, it is consistent with $\tilde {\mathcal {E}}^+\propto k_{\perp }^{-3/2}$, as would be expected at small scales once nonlinear shearing rates inevitably overwhelm reflection-related physics (see below). The evolution of $\tilde {\mathcal {E}}^{-}(k_{\perp })$ is quite different, rapidly moving to large scales at very early times. This feature is consistent with the discussion above, where we argued that the dominant scale of $\tilde {z}^{-}$ has no reason to match that of $\tilde {z}^+$, because large-amplitude $\tilde {z}^+$ eddies cause both stronger growth and stronger damping. The spectral slope rapidly reaches $\tilde {\mathcal {E}}^{-}\propto k_{\perp }^{-3/2}$ over a wide range of scales that overlaps with the range where $\tilde {\mathcal {E}}^+\propto k_{\perp }^{-1}$.

Figure 3.

Wave-action-energy spectra $\tilde {\mathcal {E}}^{\pm }(k_{\perp })$ during the imbalanced phase of the simulation. Plots (a,b) show $\tilde {\mathcal {E}}^+(k_{\perp })$ and $\tilde {\mathcal {E}}^{-}(k_{\perp })$, respectively (note the differing vertical scales), with the different colours showing different time/radii, as indicated by the colour bar. In each panel, the inset shows the best-fit power-law spectral slope, which is fit below the measured correlation scale at each $a$. Plot (c) shows both $\tilde {\mathcal {E}}^+$ (red) and $\tilde {\mathcal {E}}^{-}$ (blue) at $a\approx 5$ when the simulation is refined to the higher resolution of $n_{\perp }^{2}\times n_{z} = 8192^{2}\times 256$. Dashed black lines show various power-law slopes, highlighting a steepening of $\mathcal {\tilde {E}}^+(k_{\perp })$ at small scales (although there is not sufficient range to say whether it steepens to $\tilde {\mathcal {E}}^+\propto k^{-3/2}$ as observed in the solar wind). The inset shows the 2-D spectrum of the dominant waves $\tilde {\mathcal {E}}^+(k_{\perp },k_{z})$, illustrating how the fluctuations have decorrelated in the parallel direction (as indicated by the approximately vertical contours at larger $k_{\perp }$).

These basic features can be plausibly understood within the framework discussed above if we also consider that $\tilde {z}^{-}$ could consist of two qualitatively separate ‘classical’ and ‘anomalous’ components, as introduced in previous works Velli et al. (Reference Velli, Grappin and Mangeney1989), Verdini, Velli & Buchlin (Reference Verdini, Velli and Buchlin2009), Perez & Chandran (Reference Perez and Chandran2013) and Chandran & Perez (Reference Chandran and Perez2019). The anomalous component maintains coherence with $\tilde {z}^+$, allowing it to shear coherently over long times and thus dominating $\tilde {z}^+$'s turbulent decay. The classical part, in contrast, would be that cascaded from larger scales in $\tilde {z}^{-}$, dominating the measured spectrum but only weakly affecting the decay of $\tilde {z}^+$ because the nonlinear interactions are weak and accumulate as a random walk. Indeed, the claim above – that $\tilde {z}^{-}$ should form a $\tilde {\mathcal {E}}^{-}\propto k_{\perp }^{-3}$ spectrum due to the balance between reflections and nonlinearity – is not sustainable towards small scales. In particular, in order to form a $k_\perp ^{-3}$ spectrum, the energy injection at each scale from reflection must be larger than the flux arriving to this scale from larger scales due to nonlinear transfer. Based on the phenomenology of § 4.1 and using $\tilde {\mathcal {E}}^+\propto k_{\perp }^{-1}$, one finds that the injected flux scales as $\varepsilon ^-\propto k_\perp \tilde {z}^+(\tilde {z}^{-})^2\propto k_\perp ^{-1}$, implying that it declines towards smaller scales and will be overwhelmed by the nonlinear transfer from larger scales (Verdini et al. Reference Verdini, Velli and Buchlin2009). This idea can thus be used to motivate there being a ‘hidden’ $\tilde {\mathcal {E}}^{-}\propto k_{\perp }^{-3}$ spectrum in figure 3 that is the dominant advector of the $\tilde {z}^+$ (interestingly, the measured spectrum of the 2-D modes does follow $\tilde {\mathcal {E}}^{-}\propto k_{\perp }^{-3}$; not shown). As noted by Velli et al. (Reference Velli, Grappin and Mangeney1989) and extended to more realistic anisotropic turbulence by Perez & Chandran (Reference Perez and Chandran2013), because the $\tilde {z}^+$ cascade rate is $\varepsilon ^+\propto k_{\perp }\tilde {z}^{-} (\tilde {z}^+)^{2}$, if this is independent of $k_{\perp }$ (a constant-flux $\tilde {z}^+$ cascade), the $\tilde {z}^+$ spectrum would be $\tilde {\mathcal {E}}^+\propto k_{\perp }^{-1}$ as observed here, in previous reflection-turbulence simulations (Verdini et al. Reference Verdini, Velli and Buchlin2009; Perez & Chandran Reference Perez and Chandran2013; Chandran & Perez Reference Chandran and Perez2019; Squire, Chandran & Meyrand Reference Squire, Chandran and Meyrand2020) and in the solar wind.

4.3 Anomalous coherence

The ‘coherence assumption’ was used extensively in the discussion above in order to justify using the nonlinear time $\tau _+\sim a^{3/2}\tilde {\lambda }_{-}/\tilde {z}^{-}$ to estimate the turbulent-decay rate of the $\tilde {z}^+$ fluctuations, even though the $\tilde {z}^{-}$ fluctuations are very low amplitude and, thus, might be expected to cascade $\tilde {z}^+$ weakly. In figure 4 we diagnose this assumption numerically using the space–time Fourier spectrum (Meyrand, Galtier & Kiyani Reference Meyrand, Galtier and Kiyani2016; Lugones et al. Reference Lugones, Dmitruk, Mininni, Pouquet and Matthaeus2019), defined as

(4.8)\begin{equation} \tilde{\mathcal{E}}^{{\pm}}(k_z,\omega)=\tfrac{1}{2} \left\langle\vert \tilde{\boldsymbol{z}}^{{\pm}} (k_z,\omega)\vert^2\right\rangle_\perp,\end{equation}

where $\tilde {\boldsymbol {z}}^+(k_z,\omega )$ are the Fourier transforms in time and space of the Elsässer field. The average, $\langle \cdot \rangle _\perp$, is taken over all perpendicular wavenumbers, meaning that $\tilde {\mathcal {E}}^{\pm }(k_z,\omega )$ will be dominated by contributions from the perpendicular scales that dominate the energy spectrum at each $k_z$. In the absence of reflection, linear $\tilde {z}^{\pm }$ perturbations satisfy the dispersion relation $\omega ^{\pm }=\pm k_{z} v_{A}$, so would each show up as a single line in $\tilde {\mathcal {E}}^{\pm }(k_z,\omega )$ at $\omega = \pm k_{z} v_{A}$ (we take $k_{z}>0$). In the nonlinear simulation, the $\omega$ location of the peak of $\tilde {\mathcal {E}}^{\pm }(k_z,\omega )$ vs $k_{z}$ thus indicates the effective velocity of $\tilde {z}^{\pm }$ perturbations, while its width provides a measure of the level of nonlinear broadening due to the turbulence. Note that, because the Fourier transform is taken in $k_{z}$, rather than $k_{\|}$, care is required to ensure that the diagnostic is not affected by field-line wandering. We will see that this likely pollutes the results for $k_{z}\gtrsim 25$ in our simulations.

Figure 4.

Space–time Fourier transform (4.8) of $\tilde {\boldsymbol {z}}^+$ (a,c) and $\tilde {\boldsymbol {z}}^{-}$ (b,d). Each column is normalized to its maximum value to better illustrate the structure. Plots (a,b) show the $\chi _{\rm exp0}=960$ reflection-driven turbulence simulation at $a\approx 5$ (as in figure 2); (c,d) show the same simulation around the same time, but restarted with the reflection and expansion terms artificially removed (viz., as a normal decaying RMHD turbulence starting with initial conditions generated from the reflection-driven turbulence). While $\tilde {\boldsymbol {z}}^{-}$ fluctuations remain anomalously coherent with $\tilde {\boldsymbol {z}}^+$ in the reflection-driven simulation (a,b), the homogenous decaying turbulence does not exhibit this feature (the dominance of outwards-propagating $\tilde {z}^{-}$ fluctuations at $k_{z}\gtrsim 20$ in (d) is likely due to field-line wandering and the diagnostic should not be trusted in this range).

In figure 4(a,b) we show $\tilde {\mathcal {E}}^+(k_z,\omega )$ (a,c) and $\tilde {\mathcal {E}}^{-}(k_z,\omega )$ (b,d) in the $\chi _{A0}=96$ reflection-driven simulation. They are normalized to their maximal values at each $k_{z}$ and computed over several Alfvén crossing times around $a\approx 5$. As expected, the $\tilde {z}^+$ fluctuations concentrate in the vicinity of the Alfvén-wave predictionFootnote ⁸ $\omega \approx k_{z}v_{A}$, with modest nonlinear broadening. But, the $\tilde {z}^{-}$ fluctuations (b) are seen to propagate oppositely to linear Alfvén waves, populating the same (a,b) region as the $\tilde {z}^+$. This provides direct empirical evidence that they propagate together with $\tilde {z}^+$, leading to anomalous coherence. In the frame of the $\tilde {z}^+$ fluctuations, such $\tilde {z}^{-}$ are stationary, and can thus coherently shear the $\tilde {z}^+$ eddy over the time scale $\tau _{-}$.

To assess the role of reflection in supporting this phenomenon, in the panels (c,d) we illustrate the same plots, but for standard homogenous decaying turbulence. Specifically, we restart the reflection-driven simulation from the same time pictured in the panels (a,b), but with the reflection and expansion terms removed, then allow this turbulence to decay for several Alfvén time to measure $\tilde {\mathcal {E}}^{\pm }(k_z,\omega )$ (over this timeframe, $\tilde {z}^{-}$ decays noticeably, but $\tilde {z}^+$ does not, meaning the effect of $\tilde {z}^+$ on $\tilde {z}^{-}$ should remain similar). While $\tilde {\mathcal {E}}^+(k_z,\omega )$ (c) remains similar, we see a much wider spread in $\tilde {\mathcal {E}}^{-}(k_z,\omega )$ (d), which extends down to $\omega \approx -k_{z} v_{A}$. These general features are as expected because the $\tilde {z}^+$ modes shear the $\tilde {z}^{-}$ modes with a nonlinear time comparable to their linear time, thus forming a nonlinear frequency spread of width ${\sim }k_{z}v_{A}$. The change to $\omega >0$ dominating around $k_{z}\gtrsim 25$ is artificial, occurring because our Fourier transform in $k_{z}$ does not correctly follow the field lines, causing the measurement to be dominated by the advection of high-$k_{\perp }$ structures (presumably this same effect occurs in the left panels also, but is hidden because the fluctuations already sit at $\omega >0$).

The simplest way to understand these results is as a direct numerical demonstration of the importance of reflection in maintaining anomalous coherence in imbalanced turbulence (Chandran & Perez Reference Chandran and Perez2019). Figure 4(a,b) verifies that the $\tilde {z}^{-}$ effectively remain stationary in the frame of $\tilde {z}^+$ fluctuations; they thus do not undergo Alfvén-wave collisions and can shear $\tilde {z}^+$ coherently to enable a strong cascade. While similar ideas have appeared in a number of previous works for both homogenous and reflection-driven turbulence (Velli et al. Reference Velli, Grappin and Mangeney1989; Lithwick et al. Reference Lithwick, Goldreich and Sridhar2007; Perez & Chandran Reference Perez and Chandran2013; Chandran & Perez Reference Chandran and Perez2019; Schekochihin Reference Schekochihin2022), our results here provide a particularly clear demonstration of the effect and, via the comparison of figures 4(b) and 4(d), establish the importance of reflection in maintaining the coherence. Interestingly, Lugones et al. (Reference Lugones, Dmitruk, Mininni, Pouquet and Matthaeus2019) and Yang et al. (Reference Yang2023) have reported similar anomalous behaviour of $\tilde {z}^{-}$ in homogenous imbalanced MHD turbulence simulations with external forcing and homogeneous compressible imbalanced decaying MHD turbulence, respectively. While this does not directly disagree with our results since we have considered different physical situations (forced versus decaying and compressible versus incompressible), the topic clearly deserves more study to understand the impact of forcing (via reflection or otherwise) and compressible effects on coherence.

4.4 Wave-action anastrophy growth and the split cascade

In this section we argue that the turbulent growth of wave-action anastrophy (wave-action magnetic vector potential squared) causes $\tilde {L}_+$, the co-moving correlation scale of $\tilde {z}^+$, to rush to large scales as $\tilde {z}^+$ decays. This effect places a strong constraint on the nonlinear dynamics with interesting implications for the solar wind. It can be equivalently viewed in the expanding (physical) frame as the turbulent suppression of anastrophy decay compared with what occurs for linear waves.

4.4.1 Wave-action anastrophy

Our starting point is to note that, because $\tilde {\boldsymbol {\nabla }}_\perp \boldsymbol {\cdot }\tilde {\boldsymbol {z}}^{\pm }=0$, $\tilde {\boldsymbol {\nabla }}_\perp \boldsymbol {\cdot }\tilde {\boldsymbol {b}}_{\perp }=0$ and $\tilde {\boldsymbol {\nabla }}_\perp \boldsymbol {\cdot }\tilde {\boldsymbol {u}}_{\perp }=0$, one can define the wave-action potentials:

(4.9a–c)\begin{equation} \hat{\boldsymbol{z}}\times \tilde{\boldsymbol{\nabla}}_\perp\tilde{\zeta}^{{\pm}} = \tilde{\boldsymbol{z}}^{{\pm}},\quad \hat{\boldsymbol{z}}\times\tilde{\boldsymbol{\nabla}}_\perp\tilde{A}_z = \tilde{\boldsymbol{b}}_{{\perp}},\quad \hat{\boldsymbol{z}} \times \tilde{\boldsymbol{\nabla}}_\perp\tilde{\varPhi} = \tilde{\boldsymbol{u}}_{{\perp}}.\end{equation}

Here, $\tilde {\boldsymbol {\nabla }}_\perp$ is the co-moving-frame gradient, so these potentials differ from those naturally defined in the physical (expanding) frame, but will be more convenient here.Footnote ⁹ Equation (2.6) can then equivalently be written in terms of $\tilde {\zeta }^{\pm }$, or $\tilde {\varPhi }$ and $\tilde {A}_z$, which evolves as

(4.10)\begin{equation} \dot{a}\frac{\partial \tilde{A}_z}{\partial \ln a}+\frac{1}{a^{1/2}}\{\tilde{\varPhi}, \tilde{A}_z\} =v_{A0} \frac{\partial \tilde{\varPhi}}{\partial z} + \frac{\dot{a}}{2} \tilde{A}_z, \end{equation}

where the Poisson bracket is defined as $\{ \tilde {\varPhi },\tilde {A}_z \} = {\hat {\boldsymbol {z}}}\boldsymbol {\cdot }\tilde {\boldsymbol {\nabla }}_\perp \tilde {\varPhi }\times \tilde {\boldsymbol {\nabla }}_\perp \tilde {A}_z$ (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). Multiplying (4.10) by $\tilde {A}_z$ and integrating, we form the equation for wave-action anastrophy, $\tilde {\mathcal {A}}\equiv \langle \tilde {A}_z^2\rangle /2$:

(4.11)\begin{equation} \dot{a}\left(\frac{\partial \tilde{\mathcal{A}}}{\partial \ln a} - \tilde{\mathcal{A}}\right) = v_{A0} \left\langle \tilde{A}_z \frac{\partial \tilde{\varPhi}}{\partial z}\right\rangle = \frac{v_{A0}}{2} \left\langle \tilde{\zeta}^+ \frac{\partial \tilde{\zeta}^{-}}{\partial z}\right\rangle.\end{equation}

The nonlinear term has disappeared because anastrophy is an ideal invariant of the 2-D RMHD system, while the expansion causes $\tilde {\mathcal {A}}$ to grow (the $-\tilde {\mathcal {A}}$ on the left-hand side of (4.11)) and the three-dimensional (3-D) term $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle = - \langle \tilde {\zeta }^{-}\partial _z\tilde {\zeta }^+ \rangle$ can in principle either destroy or create it, depending on the correlation between the two Elsässer fields. Omitted in (4.11) is an additional hyper-dissipation term on its right-hand side, which can dissipate small-scale $\tilde {\mathcal {A}}$ and, thus, provide an important contribution if there exists a turbulent flux of $\tilde {\mathcal {A}}$ to small scales.

Equation (4.11) shows that if $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle$ is small in the appropriate sense, wave-action anastrophy will grow rapidly (up to $\tilde {\mathcal {A}}\propto a$), purely due to linear expansion effects.Footnote ¹⁰ As a relevant example, if the fluctuations satisfy $|\sigma _{\theta }|=1$ ($\tilde {\zeta }^+\propto \tilde {\zeta }^{-}$), lying on the edge of the circle plot in figure 1(b), then $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle =0$, driving growth of $\tilde {\mathcal {A}}$. We will now argue that in strong reflection-driven turbulence, the wave-action anastrophy grows with $a$, even in three dimensions. The argument relies on considering what occurs for propagating linear Alfvén waves, which, so long as $\varDelta =k_{z}v_{A0}/\dot {a}>1/2$ (see § 5.1), propagate with constant amplitude on average, and thus, constant $\tilde {\mathcal {A}}$. This implies that $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle$ in (4.11) must exactly balance the expansion-induced growth. Indeed, as shown in Appendix A, as an outwards ($\tilde {\zeta }^+$) fluctuation propagates, the reflected $\tilde {\zeta }^{-}$ component trails it by ${\rm \pi} /2$ in phase and has exactly the required amplitude to ensure that $v_{A0}\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle = -2\dot {a} \tilde {\mathcal {A}}$. Because the phase offset of ${\rm \pi} /2$ causes $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle$ to be as negative as possible, this implies that so long as $|\tilde {\zeta }^{-}|/|\tilde {\zeta }^+|$ remains similar to (or less than) the linear solution, any change to the phase offset between $\tilde {\zeta }^{-}$ and $\tilde {\zeta }^+$ will increase $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle$ (decrease $|\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle |$), thus causing $\tilde {\mathcal {A}}$ to grow with $a$.

For application in strong reflection-driven turbulence, it is therefore helpful to compare $\tilde {z}^{-}$ in the phenomenology of § 4.1 to what the linearly reflected $\tilde {z}^{-}$ would be for a given $\tilde {z}^+$, knowing that, if its phase offset is perfect, the latter destroys $\tilde {\mathcal {A}}$ at just the correct rate to maintain constant $\tilde {\mathcal {A}}$. The nonlinear phenomenology yields $\tilde {z}^{-}\sim \tilde {z}^+/\chi _{\rm exp}$ (§ 4.1), while the linearly reflected component is $\tilde {z}^{-} \sim \tilde {z}^+/\varDelta$ (see Appendix A). Therefore, the ratio of the two is the critical balance parameter $\chi _{A}$ – a sensible expectation given that $\chi _{A}$ is the ratio of the two effects (Alfvénic propagation and nonlinearity) that can compete with expansion to halt the growth of $\tilde {z}^{-}$. This implies that in strong ($\chi _{A}\sim 1$) reflection-driven turbulence, the amplitude of the growing $\tilde {z}^{-}$ is no larger than the amplitude needed to maintain constant $\tilde {\mathcal {A}}$. The consequence is that any modification to the linear (${\rm \pi} /2$) phase offset between $\tilde {z}^{-}$ and $\tilde {z}^+$ will decrease $v_{A0}|\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle |$ below $2\dot {a}\tilde {\mathcal {A}}$, thereby causing $\tilde {\mathcal {A}}$ to grow. While chaotic nonlinear interactions will generically act to scramble the phases of $\tilde {\zeta }^{\pm }$, we argue that reflection turbulence causes a more pronounced effect: the anomalous coherence, which leads to the high observed correlation between $-\tilde {\boldsymbol {z}}^{-}$ and $\tilde {\boldsymbol {z}}^+$ (negative $\sigma _\theta$), also precludes a large correlation between $\tilde {\zeta }^+$ and $\partial _z\tilde {\zeta }^{-}$.Footnote ¹¹ In other words, the phases are partially scrambled by the turbulence, but with a tendency for correlation between $\tilde {\zeta }^+$ and $-\tilde {\zeta }^{-}$, rather than $\tilde {\zeta }^+$ and $\partial _{z}\tilde {\zeta }^{-}$. The surprising consequence is that, while the decay rate of wave-action energy increases (up to $\tilde {E}\propto a^{-1}$) as the turbulence becomes stronger (see figure 1), the opposite is true of the wave-action anastrophy: it is approximately constant in weak turbulence (where $\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle$ remains similar to its linear value), but grows in strong turbulence.

4.4.2 The growth of $\tilde {L}_+$

From here, the arguments are standard (Kolmogorov Reference Kolmogorov1941). The wave-action energy, which is almost a true inviscid invariant during the imbalanced phase when $|\tilde {E}^{r}|\ll \tilde {E}$, decays nonlinearly due to the turbulent flux between the co-moving correlation scale $\tilde {L}_+$ and the dissipation scales. But, because the small-scale dissipation of $\tilde {\mathcal {A}}$ is proportional to the magnetic energy, for small (hyper-)viscosity, if the nonlinear dissipation of $\tilde {E}$ remains finite, the nonlinear dissipation of $\tilde {\mathcal {A}}$ must be smaller (Montgomery, Turner & Vahala Reference Montgomery, Turner and Vahala1978; Alexakis & Biferale Reference Alexakis and Biferale2018). Combined with the argument above that $|v_{A0}\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle | \lesssim 2\dot {a} \tilde {\mathcal {A}}$, we thus expect $\tilde {\mathcal {A}}$ to grow. Then, because $\tilde {E}\sim \tilde {\mathcal {A}}/\tilde {L}_+^{2}$ for imbalanced fluctuations, if $\tilde {E}$ decays while $\tilde {\mathcal {A}}$ grows (or even remains constant), this leads to remarkable phenomenon: the turbulent decay must progress with $\tilde {L}_+$ increasing rapidly in time. Specifically, taking $\tilde {\mathcal {A}}\propto a$ (assuming $|v_{A0}\langle \tilde {\zeta }^+\partial _z\tilde {\zeta }^{-} \rangle | \ll 2\dot {a} \tilde {\mathcal {A}}$ and minimal nonlinear dissipation), we find that

(4.12)\begin{equation} a^{{-}1}\tilde{E}^+ \tilde{L}_+^2 \sim \text{const.} \implies \tilde{L}_+{\propto} a, \end{equation}

where we used $\tilde {E}^+\propto a^{-1}$ from (4.4). This prediction applies to the co-moving frame, implying yet a faster increase in scales in the physical frame ($L_+\propto a^2$). Note that this law is more extreme than the standard argument for growing correlation scales in decaying 2-D MHD turbulence, which invokes only the lack of nonlinear decay of anastrophy (Hatori Reference Hatori1984; Biskamp & Welter Reference Biskamp and Welter1989). It is also worth clarifying that there is no ‘trick’ involved with the wave-action variables here: if we were instead to work in physical variables in the co-moving frame, (co-moving) anastrophy would remain constant, but the energy would linearly decay ${\propto }a^{-1}$ (and, thus, turbulently decay ${\propto }a^{-2}$) because $\boldsymbol {z}^{\pm }$ naturally decays with $a$.

The prediction (4.12) is tested in figure 5. We compute the parametric representation,

(4.13a,b)\begin{equation} X(a)={-}\dfrac{\partial \ln\tilde{z}^+_{rms}(a)}{\partial \ln a}, \quad Y(a)=\dfrac{\partial \ln\tilde{L}_+(a)}{\partial \ln a}, \end{equation}

where $\tilde {z}^+_{rms} = \sqrt {2\tilde {E}^+}$ and $\tilde {L}_+$ is computed as

(4.14)\begin{equation} \tilde{L}_{+}\equiv \int {\rm d}k_\perp\,\mathcal{\tilde{E}}^+(k_\perp)/k_\perp. \end{equation}

Here $X(a)$ and $Y(a)$ are the instantaneous scaling exponents of $1/\tilde {z}^+_{rms}$ and $\tilde {L}_+$, implying that if wave-action anastrophy, $\tilde {\mathcal {A}}\sim \tilde {E}^+\tilde {L}_+^{2}$, grows as $\tilde {\mathcal {A}}\propto a$ during the decay, then

(4.15)\begin{equation} Y(a)=X(a)+\tfrac{1}{2}.\end{equation}

This relation is independent of the decay rate of $\tilde {E}$ and, thus, the decay phenomenology. We see in figure 5 that all through the imbalanced phase ($a\lesssim 50$), $X$ and $Y$ sit almost on the line (4.15), implying $\tilde {L}_+$ grows almost as predicted by wave-action anastrophy growth (slightly more slowly). In the later dynamics, which will be described in more detail below, the fluctuation decay/growth rate ($X$) changes significantly, but wave-action anastrophy remains ${\propto }a$ as indicated by its evolution along the dotted line.

Figure 5.

Parametric representation of the instantaneous scaling exponents of $1/\tilde {z}^+_{rms}$ and the energy correlation length $\tilde {L}_+$ during the radial transport. The colours indicate the normalized radial distance $a$ (in logarithmic space). The dashed line $Y=X+1/2$ represents the theoretical expectation based on anomalous growth of anastrophy (4.12). The black star corresponds to the expected position for an anastrophy-conserving decay characterized by $\tilde {E}\propto a^{-1}$, as described in § 4.1. The black dot corresponds to the asymptotic expectation based on the linear solution (§ 5.1) for the long-wavelength expansion-dominated modes with $\varDelta <1/2$, which dominate the simulation at late times.

4.4.3 The split cascade

Physically, the fast increase in $\tilde {L}_+$ implies the energy decays through a split cascade, whereby it is forced to flow to both small and large scales simultaneously. We diagnose this surprising phenomenon directly in figure 6 by computing the Elsässer perpendicular wave-action-energy fluxes as a function of perpendicular wavenumber $k_\perp$ (Frish Reference Frisch1995):

(4.16)\begin{equation} \varPi^{{\pm}}(k_\perp)={-}a^{{-}3/2}\dfrac{2{\rm \pi}}{L_{{\perp}}}\int \dfrac{{\rm d}^3\boldsymbol{r}}{V}\left[\tilde{\boldsymbol{z}}^{{\pm}} \right]_{k_\perp}^{<}\boldsymbol{\cdot}(\tilde{\boldsymbol{z}}^{{\mp}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}_\perp\tilde{\boldsymbol{z}}^{{\pm}}), \end{equation}

where the low-pass filter is defined by

(4.17)\begin{equation} \left[\tilde{\boldsymbol{z}}^{{\pm}} \right]_{k_\perp}^{<}=\sum_{k'_{z}}\sum_{{\mid}\boldsymbol{k}'_\perp{\mid}\leqslant k_\perp}{\rm e}^{{\rm i}\boldsymbol{k}'\boldsymbol{\cdot}\boldsymbol{r}} \tilde{\boldsymbol{z}}^{{\pm}}_{\boldsymbol{k}}. \end{equation}

The split cascade of the energetically dominant field $\tilde {z}^+$ is revealed by the break between the blue and red bands that extends diagonally upwards. It is located near the measured $1/\tilde {L}_+$ at earlier times, decreasing as expected due to the conservation of anastrophy (approximately ${\propto }1/a$). On the right of the break, $\tilde {E}^+$ cascades towards small scales where reflection becomes subdominant and the hyperviscosity allows its dissipation; on the left, $\tilde {E}^+$ cascades towards large scales, allowing $\tilde {L}_+$ to increase in time. The break scale deviates modestly from the ${\propto }a^{-1}$ expectation, increasing more rapidly at early times and then slowing somewhat around $a\approx 5$ for unknown reasons, but its behaviour is broadly consistent with the evolution of $\tilde {L}_+$ (figure 5). The subdominant field $\tilde {z}^{-}$ undergoes a direct cascade during its entire evolution, aside from at the largest scales at late times, where the dynamics start becoming effectively two dimensional and balanced, differing significantly from the imbalanced phase (see below). This leads to the interesting phenomenon whereby $\tilde {z}^{-}$ and $\tilde {z}^+$ cascade in opposite directions across a modest range of intermediate scales (those above the break scale in $\varPi ^+$) during the imbalanced turbulent decay. Similar dual, counter-directional Elsässer cascades have been reported previously in flux-tube simulations of coronal holes (van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2017) and observed in high cross-helicity solar-wind streams (Smith et al. Reference Smith, Stawarz, Vasquez, Forman and MacBride2009; Coburn et al. Reference Coburn, Forman, Smith, Vasquez and Stawarz2015).

Figure 6.

Two-dimensional $k_\perp$-$a$ evolution of the Elsässer fluxes $\varPi ^+(k_\perp )$ (a) and $\varPi ^-(k_\perp )$ (b; see (4.16)). At each $a$ the $\varPi ^\pm (k_\perp )$ are normalized by their maximum over $k_\perp$ in order to better show their structure. We see clear evidence of a split cascade in $\tilde {\boldsymbol {z}}^+$, with a break between the forward and inverse cascades that migrates to larger scales with time. Although the cause of the modest deviations from the ${\propto }a^{-1}$ scaling remains unclear, the general behaviour is consistent with the discussion in the text and the evolution of the correlation length in figure 5.

5.1 Linear EBM dynamics

We will show below that by organizing itself into structures that minimize the nonlinear stresses, the system becomes effectively linear in its late stages. We thus describe basic features of the linear solution here, focusing on the difference between short-wavelength propagating (Alfvénic) waves and expansion-dominated solutions at long wavelength, which grow continuously with $a$. These linear solutions are illustrated in figure 7, starting from pure $\tilde {\boldsymbol {z}}^+$ fluctuations in the initial conditions. Their characteristics, including the growth of expansion-dominated modes, have been studied using various methods in global geometries in a number of previous works (Heinemann & Olbert Reference Heinemann and Olbert1980; Hollweg Reference Hollweg1990; Hollweg & Isenberg Reference Hollweg and Isenberg2007); they are not an artefact of the EBM.

Figure 7.

Solutions of the linearised equations (5.1), starting from the initial condition $\tilde {z}^{-}(0)=0$ and $\tilde {z}^+(0)=\sqrt {2}$ with different values of $\varDelta$ as labelled. Solid lines show $|\tilde {z}^+(a)|$; dotted lines show $|\tilde {z}^{-}(a)|$. Here $\varDelta >1/2$ modes (red, yellow and green curves), which are dominated by Alfvénic forces, exhibit wave-like behaviour with no long-term growth or decay (${\rm Im}(\omega )=0$); $\tilde {z}^+$ propagates Alfvénically with an oscillating phase and approximately constant amplitude, while the amplitude of $\tilde {z}^{-}$ alternates up and down over the wave period as $\tilde {z}^{-}$ moves in and out of phase with the reflection forcing from $\tilde {z}^+$ (its maximum amplitude scales ${\propto }\varDelta ^{-1}$; see (5.3) and Appendix A). In contrast, long-wavelength $\varDelta <1/2$ modes with ${\rm Re}(\omega )=0$ (blue and black curves) do not oscillate like waves at all because the reflection overwhelms the Alfvénic restoring force (see (5.4); Heinemann & Olbert Reference Heinemann and Olbert1980). The amplitude of the magnetically dominated mode grows as $|\tilde {z}^{\pm }(a)|\propto a^{|\omega ^\pm |}$, with the growth rate $|\omega ^\pm |=\tfrac 12\sqrt {1-4\varDelta ^2}$ depending only weakly on $\varDelta$ (cf. blue and black curves).

The full linear solution is easily obtained by ignoring the nonlinear terms in (2.7) and assuming divergence-free plane-wave solutions of the form $\tilde {\boldsymbol {z}}^{\pm } = \tilde {z}^{\pm }(a) \exp ({{\rm i}k_\perp y + {\rm i}k_z z}) \hat {\boldsymbol {x}}$. This gives

(5.1)\begin{equation} \dfrac{\partial \tilde{z}^{{\pm}}}{\partial \ln{a}}+\begin{pmatrix} {\rm i}\varDelta & 1/2 \\ 1/2 & -{\rm i}\varDelta \end{pmatrix}\begin{pmatrix} \tilde{z}^+\\\tilde{z}^{-}\end{pmatrix}=0,\end{equation}

where $\varDelta = k_{z}v_{A}/(\dot {a}/a) = k_{z}v_{A0}/\dot {a}$ (using the time variable $\ln a$ eliminates explicit time dependence from the linear system), allowing one to insert the ansatz,

(5.2)\begin{equation} \tilde{z}^{{\pm}}(\ln a) = \tilde{z}^{{\pm}}_w \exp({\rm i}\omega \ln a), \end{equation}

where $\tilde {z}^{\pm }_w$ is the complex amplitude of $\tilde {z}^{\pm }(a)$. The general solution to (5.1) can then be formed via the eigenmodes,

(5.3)\begin{equation} \xi^{{\pm}}=\tfrac{1}{2}\tilde{z}^{{\pm}}_w \pm i\left(\varDelta\mp\omega^{{\pm}}\right)\tilde{z}^{{\mp}}_w, \end{equation}

which evolve as $\xi ^\pm (a) = \xi ^\pm _0\exp ({\rm i}\omega ^\pm \ln a)$ from initial conditions $\xi ^\pm _0$, where the eigenfrequencies $\omega ^\pm$ are

(5.4)\begin{equation} \omega^{{\pm}}={\pm} \sqrt{\varDelta^2-1/4}. \end{equation}

We see that $\varDelta =1/2$ marks the boundary between oscillating Alfvénic modes and growing (or decaying) expansion-dominated modes: for $\varDelta >1/2$, $\omega ^{\pm }$ is real and $\tilde {z}^{\pm }$ oscillates with frequency $\omega ^\pm$, albeit with a minority-reflected $\tilde {z}^{\mp }$ component that inevitably accompanies any $\tilde {z}^{\pm }$ fluctuation; for $\varDelta <1/2$, $\omega ^{\pm }$ is imaginary and modes grow exponentially, $\tilde {z}^{\pm } \propto \exp ({|\omega ^\pm |\ln a}) = a^{|\omega ^\pm |} = a^{\sqrt {1-4\varDelta ^2}/2}$, because the expansion overwhelms the Alfvénic restoring force. The growing expansion-dominated mode, with $\omega = i\sqrt {1/4-\varDelta ^{2}}$, is magnetically dominated with $\tilde {\boldsymbol {z}}^{-}\approx -\tilde {\boldsymbol {z}}^+$ and $|\tilde {\boldsymbol {b}}_{\perp }|\gg |\tilde {\boldsymbol {u}}_{\perp }|$, while the decaying mode (${\rm Im}(\omega )<0$) is $\tilde {\boldsymbol {u}}_{\perp }$ dominated. Physically, the ${\sim }a^{1/2}$ growth of $\tilde {z}^{\pm }$ corresponds to $|\boldsymbol {B}_{\perp }|\propto a^{-1}$ ($|\boldsymbol {B}_{\perp }|/|\bar {\boldsymbol {B}}|\propto a$) so that $\boldsymbol {b}_{\perp }=\boldsymbol {B}_{\perp }/\sqrt {4{\rm \pi} \rho }$ is constant (Velli, Grappin & Mangeney Reference Velli, Grappin and Mangeney1991; Grappin et al. Reference Grappin, Velli and Mangeney1993). Clearly, if there exists any power in such expansion-dominated modes at early times, they will inevitably come to dominate the late-time evolution, overtaking the Alfvénic ($\varDelta >1/2$) modes.

In our simulations with $\varDelta _{\rm box}=10$ (2.12), only the $k_{z}=0$ periodic mode lies in this expansion-dominated regime. But, the properties of expansion-dominated modes are rather insensitive to $k_{z}$ for $\varDelta <1/2$: the modes have no real frequency (oscillating) part and growth rates that exhibit only a small correction compared with the $\varDelta =0$ mode (${\rm Im}(\omega ^+)\approx 1/2-\varDelta ^{2}$ for small $\varDelta$). Therefore, we argue that their dynamics should be adequately captured by the simulation, even though true $k_{z}=0$ modes are obviously not possible in a realistic non-periodic system. In reality, if we assume that the longest-wavelength modes possible are those of the system scale, with $k_{z}\sim 1/R$, then the minimum $\varDelta$ available to the system is $\varDelta _{\rm min}\sim (v_{A}/R)/(U/R)\sim v_{A}/U<1$. Thus, in the super-Alfvénic ($v_{A}< U$) wind it is always consistent to assume that the expansion-dominated modes exist, and indeed, the range of such modes available to the system will be an increasing function of radius. This feature, whereby $\varDelta _{\rm min}$ decreases with radius, is clearly not possible to capture in the EBM with a fixed $L_{z}$ (it is captured by global linear solutions; Heinemann & Olbert Reference Heinemann and Olbert1980; Hollweg & Isenberg Reference Hollweg and Isenberg2007), so the impact of this physics should be tested in future work using flux-tube simulations.

5.2 Transition into the magnetically dominated phase

During the imbalanced phase in our simulations, the turbulence appears to remain strong with $\chi _{A}\sim 1$, rapidly adjusting its parallel correlation length $\ell _{\|}$ towards larger scales as the turbulence decays (after its initial transient adjustment from the initial conditions, which occurs by $a\approx 1.2$). This phase ends and the decay deviates from the $\tilde {E}\propto a^{-1}$ phenomenology of § 4.1, once it decays sufficiently so that the box-wavelength modes ($k_{z}=2{\rm \pi} /L_{z}$) become weak ($\chi _{A}<1$). This occurs around $a\approx 25$ for the solid lines in figure 1, which agrees well with the value expected from solving $z^+/\tilde {\lambda }_+\simeq v_{A} 2{\rm \pi} /L_{z}$ with $\tilde {z}^+\propto a^{-1/2}$ and $\tilde {\lambda }_+\propto a^{-1/2}$. Following this, the expansion-dominated modes inevitably take over, driving the system towards the $|\tilde {\boldsymbol {b}}_{\perp }|\gg |\tilde {\boldsymbol {u}}_{\perp }|$ linear solution that grows with $|\tilde {\boldsymbol {b}}_{\perp }|\propto a^{1/2}$.

More generally, without the limitations of our periodic box, this transition should be understood by noting that if the turbulence remains strong with $\chi _{A}\sim 1$ throughout its decaying imbalanced phase (as appears to be the case until it becomes artificially constrained by the box), the transition to the balanced regime, at $\chi _{\rm exp}\sim 1$, will occur when $\varDelta =\chi _{\rm exp}/\chi _{A}\sim 1$, viz., at the same time that the dominant modes in the system become expansion dominated. This pleasing consistency of the phenomenology argues that the system cannot reach the balanced phase while still dominated by Alfvénic physics and suggests that the large scales in the balanced phase will not be critically balanced in the usual sense (because their linear physics is dominated by expansion not Alfvénic propagation). This property should hold so long as the turbulence remains strong during the imbalanced decay phase and transitions into the balanced phase at $\chi _{\rm exp}\sim 1$ (i.e. independently of the evolution of $\tilde {\lambda }_+$ or other uncertainties in § 4.1).

Following this transition, any further turbulent decay will tend to increase the nonlinear time, thus driving the system inevitably towards the linear regime where expansion dominates both Alfvénic and nonlinear effects. This can be seen by noting that unless $\tilde {\lambda }_+$ decreases, then even the fastest-possible linear growth, $\tilde {z}^+\sim \tilde {z}^{-}\propto a^{1/2}$, leads to $\chi _{\rm exp}\sim a^{-1/2} (\tilde {z}^+/\tilde {\lambda }_+)/\dot {a}$ remaining constant (the dominance of expansion over Alfvénic propagation is guaranteed because $\varDelta \lesssim 1$). However, we see from the perpendicular structure shown in figure 2 that the system approaches this expansion-dominated state in an interesting and non-trivial way: rather than simply decaying to low amplitudes to reduce the nonlinear time, it organizes itself into isolated, coherent structures that approach nonlinear solutions in which the magnetic tension balances the pressure. This self-organization thus defeats prematurely the nonlinear couplings and turbulent dissipation, precipitating the system into magnetically dominated Alfvén vortices that behave almost linearly.

Because the system becomes expansion dominated with little turbulent dissipation, its growth must also satisfy the prediction of (4.15) for the growth of wave-action ansatrophy $\tilde {\mathcal {A}}\propto a$. Thus, during the transition as it moves into the balanced phase, the system evolves downwards along the $Y=X+1/2$ line in figure 5, tending towards the point $X=-1/2$, $Y=0$ that characterizes purely linear evolution.

5.3 Emergence of Alfvén vortices

At this point, the story is mostly finished as far as the turbulent heating and dissipation is concerned: as the system becomes balanced, it also starts shutting off its nonlinear dissipation, creating long-parallel-wavelength perpendicular structures that grow with $|\tilde {\boldsymbol {b}}_{\perp }|\sim a^{1/2}$. However, the quasi-circular structures that emerge (see figure 2) are of significant interest, both for comparison to in-situ observations, and because they are picturesque illustrations of the ‘cellularization’ of turbulence (Matthaeus et al. Reference Matthaeus, Wan, Servidio, Greco, Osman, Oughton and Dmitruk2015) – a vivid example of self-organization (Horiuchi & Sato Reference Horiuchi and Sato1985). They can be understood using a classical variational argument (Matthaeus & Montgomery Reference Matthaeus and Montgomery1980). Motivated by the turbulent wave-action anastrophy growth, we minimize the Alfvénic magnetic energy per unit volume, $\langle |\boldsymbol {b}_\perp |^2\rangle /2=a^{-1}\tilde {E}^b,$ subject to the constancy of the anastrophy per unit volume, $\langle A_z^2\rangle /2=a^{-1}\tilde {\mathcal {A}}$ (by this choice of variables, we factor out the expansion-induced dependence on $a$; both $\langle |\boldsymbol {b}_\perp |^2\rangle$ and $\langle A_z^2\rangle$ remain constant under linear evolution for $a\gg 1$). Such minimization requires that during the relaxation process the kinetic energy is dissipated completely, leaving a pure magnetic state. It is thus aided significantly by expansion, which damps $\boldsymbol {u}_{\perp }$ but not $\boldsymbol {b}_{\perp }$ (see (2.3) and (2.2)) and preferentially increases the energy content of the longest-parallel-wavelength modes (thus creating quasi-2-D dynamics at large radial distances).

These arguments lead us to the variational problem

(5.5)\begin{equation} a^{{-}1}\delta \int {\rm d}^{3}\boldsymbol{r}\left( |\tilde{\boldsymbol{\nabla}}_\perp \tilde{A}_z|^2-\varLambda \tilde{A}_z^2\right)=0, \end{equation}

where $\delta$ denotes the functional derivative and $\varLambda$ a Lagrangian multiplier. Identifying $\varLambda$ with a characteristic scale $K_\perp$ via $\varLambda =-K_{\perp }^2$, the Euler–Lagrange equation becomes the Helmholtz equation,

(5.6)\begin{equation} \tilde{\nabla}^{2}\tilde{A}_z={-}K_{{\perp}}^2 \tilde{A}_z. \end{equation}

Recalling that $\tilde {A}_z$ evolves as a passive scalar in two dimensions (see (4.10)), we now imagine some region, or ‘cell’, in the domain that can change shape and mix in order to approach the minimum energy state, viz., the solution of (5.6) with the minimum possible $K_\perp$. The argument is effectively that the Lagrange multiplier $K_\perp$ should be piecewise constant, enforcing the minimization principle across patch-like cells where the turbulence becomes suppressed. We assume the value of $\tilde {A}_z$ outside the cell in question to be approximately constant (based on figure 2, this may not be as unreasonable as it sounds), which fixes some boundary condition $\tilde {A}_z = A_B$ on its edge. The area of the cell must remain constant because the $\tilde {\boldsymbol {u}}_{\perp }$ that advects $\tilde {A}_z$ is incompressible (similarly, the average of $\tilde {A}_z$ across the cell is fixed) – we are therefore interested in a solution of (5.6) that is as compact as possible for a given $K_\perp$, thereby providing the lowest energy (smallest $K_\perp$) for a given sized cell. This is afforded by a cylindrically symmetric cell, so we define $(r,\theta )$ as the polar coordinates centred on the cell, yielding $\tilde {A}_z\propto J_{0}(K_\perp r)$ as the only $\theta$-independent solution that does not diverge as $r\rightarrow 0$. Note that an arbitrary constant can be added to the solution by changing the gauge of $\tilde {A}_z$, but this must be added directly into the original variational problem.

Collecting these constraints, we obtain the perfectly circular magnetic-vortex solution

(5.7)\begin{equation} \begin{cases} \tilde{A}_z(r)=A_{0}J_{0}(K_\perp r), & r< r_c, \\ \tilde{A}_z(r)= A_B, & r\geqslant r_c, \end{cases} \end{equation}

where $r_c$ is the radius of the cell, at which $A_{0}J_{0}(K_\perp r)=A_B$ (as required to satisfy the boundary condition). Note that the two constants $A_0$ and $K_\perp$ are determined through the fixed area of the cell, the initial wave-action anastrophy and the boundary conditions (assuming $\tilde {A}_z$ is continuous at the start of the relaxation this will not provide a third constraint). This leaves no freedom to allow the first derivative of $\tilde {A}_z$ (i.e. the magnetic field) to be continuous, leading to an inevitable $\tilde {\boldsymbol {b}}_{\perp }$ discontinuity across the cell boundary and a strong ring of current surrounding the cell. These features are clearly observed in the simulation, as shown in figure 8, where we zoom in on various observed cells and highlight the large boundary currents (c).

Figure 8.

(a) Snapshot of the magnetic field modulus in a plane perpendicular to $B_{0}$ at $a=250$. (b) Close-up corresponding to the marked region on the left, illustrating Alfvén vortices colliding and merging through reconnection. The black circles mark the regions over which azimuthal averages have been computed to fit the Alfvén vortex solution (5.9) in figure 9. (c) Same region as the (b), but showing the out-of-plane current. This reveals sets of intense current rings, a hallmark of the ground state Alfvén vortices.

The solution (5.7) corresponds to a particular case of so-called Alfvén vortex solutions, in particular the vortex ‘monopole’ (Petviashvili & Pokhotelov Reference Petviashvili and Pokhotelov1992; Alexandrova Reference Alexandrova2008; Lion, Alexandrova & Zaslavsky Reference Lion, Alexandrova and Zaslavsky2016). As well as resulting from the variational argument, they arise as nonlinear solutions of ideal incompressible MHD equations. Indeed, the Helmholtz equation (5.6) can instead be obtained by assuming zero flow $\boldsymbol {u}_{\perp }=0$ and $k_z\ll k_\perp$, which gives, from the momentum equation (2.6),

(5.8)\begin{equation} \lbrace \tilde{A}_z,\tilde{\boldsymbol{\nabla}}_\perp^2 \tilde{A}_z \rbrace = 0. \end{equation}

Any functional relation $\tilde {\boldsymbol {\nabla }}_\perp ^2 A_z = f(A_z)$ satisfies (5.8), which subsumes any solution of the Helmholtz equation (5.6) and, thus, (5.7). In this minimum energy, constant-anastrophy solution, the contours of $A_z$ and $\tilde {\boldsymbol {\nabla }}_\perp ^2 A_z$ are circularly symmetric with aligned gradients, thus nullifying the Poisson bracket (5.8) (Grošelj et al. Reference Grošelj, Chen, Mallet, Samtaney, Schneider and Jenko2019). This nonlinear solution involves the magnetic tension balancing the perpendicular pressure.

The theoretical considerations presented above provide more than a qualitative explanation for the turbulent cellularization observed. We fit the magnetic eddies highlighted in figure 8 using the functional form

(5.9)\begin{equation} \tilde{A}_z(r)=\tilde{A}_{z0}J_{0}(K_\perp r)\left[1-f(r)\right] +\tilde{A}_z(r_c), \end{equation}

where $f(r)$ is the logistic function $f(r)=(1+\exp ({-\kappa (r-r_c)}))^{-1},$ which is effectively a step function that accounts for finite diffusive effects through the ‘logistic steepness’ parameter $\kappa$. The result of such a fit is shown in figure 9 and demonstrates that the structures observed are unequivocally the minimum-anastrophy Alfvén vortices (5.7).

Figure 9.

(a) Comparison between the absolute value of the magnetic vector potential obtained from numerical simulation at $a=250$ (coloured lines) and the analytical prediction (5.9) (black dots). The red, blue and green lines have been obtained from the Alfvén vortices labelled 1, 2 and 3 after an azimuthal average about their centre (denoted $|\langle \tilde {A}_z\rangle _\theta |$). The inset represents the analytical prediction for the magnetic field, highlighting the presence of a discontinuity at the vortex boundary. (b,c) Two-dimensional representation of the solution (5.9) for the magnetic field modulus $| \boldsymbol {b}_\perp |$ and the absolute value of the magnetic current $| j_z|$ (the colour scales are arbitrary).

In figure 10 we illustrate the magnetic- and kinetic-energy spectra, $\tilde {\mathcal {E}}^{b}$ and $\tilde {\mathcal {E}}^{u}$, respectively. The strong magnetic dominance at large scales leads to a steeper magnetic spectrum, approximately $\tilde {\mathcal {E}}^{b}\propto k_{\perp }^{-5/3}$ at large scales, with a flatter velocity spectrum that eventually joins the magnetic spectrum at small scales. This is qualitatively similar to those observed at large scales during very low cross-helicity periods in the solar wind (Tu & Marsch Reference Tu and Marsch1991). The inset shows the 2-D $k_{\perp },k_{z}$ magnetic-energy spectrum, illustrating the dominance of the $k_{z}=0$ 2-D modes at these late times. The velocity fluctuations seem to be dominated by regions between the individual magnetic cells, arising from the coalescence of the Alfvén vortices through magnetic reconnection, which generate outflows in the reconnection exhausts. The Alfvén vortices thus slowly move around, thereby generating further collisions. As the simulation progresses, ever larger magnetic structures are generated via mergers of Alfvén vortices, creating further outflows that trigger more merging, thus minimizing the total energy and causing a slow nonlinear decay (this is overwhelmed by the expansion-induced growth). However, this hierarchical process, which is the basis of 2-D MHD decaying turbulence theories (Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019), is impeded by the expansion, which acts as an additional damping of the outflows, hindering the nonlinear dissipation; indeed, the turbulent dissipation rate in our simulations, which is measured to scale as ${\sim }a^{-0.2}$ at late times, is slower than in 2-D MHD (Hatori Reference Hatori1984; Zhou et al. Reference Zhou, Bhat, Loureiro and Uzdensky2019). At large radial distances, these nonlinear effects therefore tend to ‘freeze up’ and the structures become more and more static in time, growing at almost the rate predicted from linear theory (${\propto }a^{1}$).

Figure 10.

Wave-action magnetic-energy spectrum $\tilde {\mathcal {E}}^{b}$ and kinetic-energy spectrum $\tilde {\mathcal {E}}^{u}$ at $a=250$ (cf. figure 2e,f). The magnetic energy significantly dominates at large scales, with a steeper slope that eventually joins the velocity spectrum at small scales. The inset shows the 2-D $k_{\perp },k_{z}$ spectrum of magnetic energy, illustrating how it is significantly dominated by the 2-D $\varDelta =0$ modes (the only expansion-dominated $\varDelta <1/2$ modes in our domain).

The stability of such structures is an interesting question that we do not study in detail. The inevitability of the intense current rings suggests that at a sufficiently high Lundquist number these ‘ground state’ Alfvén vortices will become tearing-unstable and break up into plasmoid chains confined on rotating rings. Allowing for the finite length of the structures and/or compressible effects may also lead to instabilities. Indeed, given the background mean field, the equilibrium (5.7) is effectively a screw pinch, with its nonlinear equilibrium resulting from the balance between the curvature/tension force of $\tilde {\boldsymbol {b}}_{\perp }$ and the pressure gradient. Such equilibria can be unstable to kink instability and sausage instabilities (Kruskal, Schwarzschild & Chandrasekhar Reference Kruskal, Schwarzschild and Chandrasekhar1954; Schuurman, Bobeldijk & de Vries Reference Schuurman, Bobeldijk and de Vries1969). Thus, by assuming the plasma to be incompressible with constant density, the RMHD model may miss important effects in their description, particularly instabilities that could aid in their destruction and enhance nonlinear dissipation.

6.1 Two-scale transport theory

Two-scale transport theory splits fields into a large-scale, quasi-stationary compressible part and a small-scale incompressible, fluctuating component. This method introduces a classical closure problem that requires a set of approximations to truncate the nonlinear terms arising from triple correlations. For all works discussed below, one-point closures are used, which, together with assumptions of particular turbulence symmetries, simplify the system sufficiently to allow the derivation of a tractable theory. The landmark studies of this approach are Marsch & Tu (Reference Marsch and Tu1989), Zhou & Matthaeus (Reference Zhou and Matthaeus1989), Velli et al. (Reference Velli, Grappin and Mangeney1989) followed by Matthaeus et al. (Reference Matthaeus, Oughton, Pontius and Zhou1994) and Matthaeus, Zank & Oughton (Reference Matthaeus, Zank and Oughton1996). All of these works include the effect of spherical expansion and, therefore, reflection, but there exists a divide between models that focus on inter-stream shear and compressions as the primary driver of turbulence (termed ‘shear-driven models’ below), and those – like this work – that focus on Alfvénic physics (termed ‘reflection-driven models’ below). There is, of course, an overlap between the two, particularly for more sophisticated models that include a wider range of physical effects (e.g. van der Holst et al. Reference van der Holst, Sokolov, Meng, Jin, Manchester, Tóth and Gombosi2014; Usmanov et al. Reference Usmanov, Matthaeus, Goldstein and Chhiber2018).

6.2 Numerical experiments

In conjunction with the development of turbulent transport models, direct simulations of turbulence are crucial for testing closure assumptions and revealing unexpected dynamical effects. Two different approaches have been considered: the Eulerian (flux tube) and the Lagrangian (expanding box). The Eulerian approach considers turbulence in a fixed narrow magnetic flux tube centred on a radial magnetic field line extending outward from the Sun (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Matsumoto & Suzuki Reference Matsumoto and Suzuki2012; Perez & Chandran Reference Perez and Chandran2013; Shoda et al. Reference Shoda, Suzuki, Asgari-Targhi and Yokoyama2019). The Lagrangian approach considers the evolution of a plasma volume advected by a spherical wind with constant speed, the so-called EBM, as adopted in this paper (Grappin et al. Reference Grappin, Velli and Mangeney1993; Dong et al. Reference Dong, Verdini and Grappin2014). The two approaches are closely related. In super-Alfvénic regions, in the limit in which the weaker field $\boldsymbol {z}^{-}$ is dissipated before it can propagate the length of the box (that is, for $\chi \geqslant 1$), the Eulerian and Lagrangian approaches should give similar results (aside from subtleties related to the periodic domain (see § 5.1). This implies that our results should be similar to the large-radius limit of RMHD flux-tube wind models such as Perez & Chandran (Reference Perez and Chandran2013) and van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2016). However, as far as we are aware, all such works have focused on the sub-Alfvénic regime, and so it is not possible to make any detailed comparison at this time. This should be a priority for future work.

Our results are more directly comparable to previous expanding-box literature, although there exist important differences in the model, initial conditions and analysis. As concerns the model, a clear difference is that most previous EBM works have used compressible MHD; by using RMHD, we have neglected compressibility and large-amplitude effects (Verdini et al. (Reference Verdini, Velli and Buchlin2009, Reference Verdini, Grappin, Pinto and Velli2012) also use simplified models). Understanding how these might influence the results is an important area for future work that requires wrestling with reflection-driven turbulence of spherically polarized fields and switchbacks (Squire et al. Reference Squire, Chandran and Meyrand2020; Johnston et al. Reference Johnston, Squire, Mallet and Meyrand2022). Another important difference is that most previous works (e.g. Dong et al. Reference Dong, Verdini and Grappin2014; Verdini & Grappin Reference Verdini and Grappin2015, Reference Verdini and Grappin2016; Montagud-Camps et al. Reference Montagud-Camps, Grappin and Verdini2018) have used initial states with $\sigma _r\sim 0$ and $\sigma _c\sim 0$, whereas we initialize with Alfvénic states, implying that the $\tilde {z}^{-}$ fluctuations are all created by reflections in our simulations. We believe this fully imbalanced initial condition to be the most relevant choice, since $\sigma _r,\sigma _c\sim 0$ turbulence is not generated by reflection and is not commonly observed in the solar wind either close to the Sun, where $\sigma _c$ is almost always observed to be ${\simeq }1$, or at larger radii (Bruno & Carbone Reference Bruno and Carbone2013; Wicks et al. Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013). The difference in initial conditions is important to the dynamics, since it leads to a turbulent decay that remains confined near the edge of the $\sigma _r^2+\sigma _r^2=1$ circle throughout the radial transport, enabling us to identify the two distinct regimes that characterize the evolution of plasma turbulence in a radially expanding flow. Montagud-Camps et al. (Reference Montagud-Camps, Grappin and Verdini2020) and Grappin et al. (Reference Grappin, Verdini and Müller2022) have also studied the expanding turbulence evolution from large $\sigma _c$, finding, similar to us, that the cross-helicity decreases while creating negative residual energy as the turbulent cascade develops (see also Squire et al. Reference Squire, Chandran and Meyrand2020). Consistent with our results, they also find weaker turbulence at lower $\chi _{A}$ and slower decay at lower $\chi _{\rm exp}$, and curiously, that the decay rate of balanced turbulence is similar to imbalanced turbulence when $\chi _{\rm exp}> 1$. A final significant difference with all previous EBM works is that we have considered initial conditions for which the correlation length is far from the box scale; this allows the development of the split cascade, the growth of the integral scale and the formation of the $1/f$ range, a key part of our results. In contrast, most simulations in Squire et al. (Reference Squire, Chandran and Meyrand2020), Montagud-Camps et al. (Reference Montagud-Camps, Grappin and Verdini2020) and Grappin et al. (Reference Grappin, Verdini and Müller2022), with imbalanced initial conditions near the box scale, exhibit spectral scalings $\propto k_\perp ^{-3/2}$ as in standard MHD turbulence. It should be a priority of future work to understand these scalings more generally.

7.1 Parameters and scales of the solar wind

Let us first remind the reader that our simulation parameters were chosen primarily to test and understand the dynamics of reflection-driven turbulence, rather than simulate specific solar-wind streams. In particular, the large $\chi _{\rm exp0}$ and small-scale $\tilde {z}^+$ (small $\tilde {L}_+/L_\perp$) are more extreme than occurs in the super-Alfvénic solar wind, while the distance range (up to $a=1000$) is wider than regularly considered in observational studies. Via the phenomenology of § 4.1, this leads to a longer phase of imbalanced-turbulence decay before the transition to the balanced phase, thus improving our analysis of these dynamics. Moreover, our highly idealized initial conditions are clearly inappropriate, since significant evolution will have occurred before reaching $R_{A}$ in the sub-Alfvénic regions of the wind.

We can estimate more realistic parameters using recent Parker Solar Probe results from $R\simeq R_{A}$ (Kasper et al. Reference Kasper2021). Using $L_+\sim U/f_{\rm out}$, where $f_{\rm out}$ is a characteristic measured frequency at the energy-dominant scale (we take the spectral break in figure 2 of Kasper et al. Reference Kasper2021), one finds that

(7.1)\begin{align} \chi_{\exp} & = \frac{z^+_{rms}/\tilde{L}_+}{\dot{a}/a} \sim f_{\rm out} \frac{R}{U^2} {z^+_{rms}} \nonumber\\ & \simeq 63 \,\mathcal{M}_{A}^{{-}1}\frac{f_{\rm out}}{2\times 10^{{-}3}\ {\rm Hz}} \frac{R/(18 R_\odot)}{U/400\ {\rm kms}^{{-}1}}\frac{z^+_{rms}}{v_{A}}. \end{align}

This estimate ignores many uncertainties, including the difference between parallel and perpendicular scales, differences between streams and the violation of Taylor's hypothesis near the Alfvén point, but should at least give an order-of-magnitude estimate of the $\chi _{\rm exp}$ of the $z^+$ fluctuations that dominate the total energy. We see that for observed fluctuation amplitudes, we expect $\chi _{\rm exp}\gg 1$, but not nearly so large as the $\chi _{\rm exp0}$ chosen for our most extreme simulation. This further implies that the transition into the balanced, magnetically dominated regime will occur at smaller radii than implied by figure 1, and the separations between the scales of $z^+$ and $z^{-}$, or between the initial and final scales of $z^+$, will be much smaller.

7.2 Relation to specific observations in the solar wind

Here we outline various predictions of reflection-driven turbulence that can be directly compared with solar-wind observations. We particularly focus on the importance of $\chi _{\rm exp}$, including the possibility that a correlation of $\chi _{\rm exp}$ with wind speed could naturally explain numerous other well-known observational correlations. While we reference various previous works throughout this discussion, a more detailed description of how this work fits into previous literature and ideas is given above in § 6.

7.2.1 Imbalance evolution

There has been substantial previous literature devoted to understanding the observed evolution of imbalance (normalized cross-helicity) with radius, as well as its correlation with wind speed (Roberts et al. Reference Roberts, Goldstein, Klein and Matthaeus1987; Tu & Marsch Reference Tu and Marsch1995). A particular focus has been understanding why the imbalance is observed to decrease with increasing radius in the solar wind, even though decaying MHD turbulence simulations (and theory) robustly show (and predict) the opposite (Dobrowolny, Mangeney & Veltri Reference Dobrowolny, Mangeney and Veltri1980; Grappin et al. Reference Grappin, Frisch, Pouquet and Leorat1982; Chen et al. Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011; Schekochihin Reference Schekochihin2022). Some suggestions invoke interactions between different streams as the dominant influence (Roberts et al. Reference Roberts, Goldstein, Matthaeus and Ghosh1992; Matthaeus et al. Reference Matthaeus, Minnie, Breech, Parhi, Bieber and Oughton2004; Adhikari et al. Reference Adhikari, Zank, Bruno, Telloni, Hunana, Dosch, Marino and Hu2015), or parametric decay of Alfvén waves (Goldstein Reference Goldstein1978), but we see that the imbalance decrease is naturally explained by reflection without invoking any other physics. While we are certainly not the first to suggest this (e.g. Velli et al. Reference Velli, Grappin and Mangeney1989; Zhou & Matthaeus Reference Zhou and Matthaeus1990; Oughton & Matthaeus Reference Oughton and Matthaeus1995; Verdini et al. Reference Verdini, Velli and Buchlin2009; Grappin et al. Reference Grappin, Verdini and Müller2022), our simulations and phenomenology clarify why this occurs and provide simple, testable predictions that (to the best of our knowledge) have not appeared in previous literature.

As argued above, the key parameter governing the imbalanced decay phase is $\chi _{\rm exp}$, the ratio of nonlinear to expansion rates. The basic phenomenology of § 4.1 (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Verdini & Velli Reference Verdini and Velli2007) predicts $\tilde {z}^{-} \sim \tilde {z}^+/\chi _{\rm exp}$, with $\chi _{\rm exp}\sim (z^+/\lambda _+)/(\dot {a}/a)$ seen to scale as $\chi _{\rm exp}\propto a^{-3/2}$ in our simulations where $\tilde {E}^+\sim a^{-1}$ and $\tilde {\lambda }_+\sim a^{1/2}$ (as inferred from the evolution of $\tilde {z}^{-}$). This implies that $\sigma _{c}$ evolves as

(7.2)\begin{equation} \sigma_{c} \sim \begin{cases} \dfrac{1-\chi_{\rm exp}^{{-}2}}{1+\chi_{\rm exp}^{{-}2}}, & \chi_{\rm exp}\gtrsim 1, \\[10pt] 0 , & \chi_{\rm exp},\lesssim 1, \end{cases} \end{equation}

which, for $\chi _{\rm exp} \sim \chi _{\rm exp0}a^{-3/2}$, stays in a highly imbalanced state near $\sigma _{c}= 1$ across a wide range of $a$ before rapidly dropping towards zero around the radius $a\sim \chi _{\rm exp0}^{2/3}$ (the exact power-law exponent $-3/2$ makes no difference to this basic picture). The model thus naturally explains the observed radial dependence of the turbulence from imbalanced to balanced. Similarly, the observed differences between fast and slow streams would be well explained if fast-wind streams start with larger $\chi _{\rm exp0}$ around $R_{A}$ (and, therefore, also maintain larger $\chi _{\rm exp}$ throughout their evolution). This prediction is easily testable observationally, aside from potential ambiguity of the outer-scale definition in $\chi _{\rm exp}$. It is also physically expected based on reflection-driven models of the sub-Alfvénic regions (Cranmer et al. Reference Cranmer, van Ballegooijen and Edgar2007; Chandran Reference Chandran2021), in which slower streams arise because more of the outward-fluctuation energy is dissipated at low altitudes (Halekas et al. Reference Halekas2023), thus giving a lower amplitude at large radii and, therefore, a lower $\chi _{\rm exp}$. For the stream with $\chi _{\rm exp0}\simeq 60$ discussed above ((7.1); Kasper et al. Reference Kasper2021), evolution with $\chi _{\rm exp}\propto a^{-3/2}$ predicts that the fluctuations should reach small $\sigma _{c}$ around $a\approx 15$, or a little beyond $1\ {\rm AU}$ – certainly not unreasonable.

7.2.2 The $\sigma _c-\sigma _r$ ‘circle plot’

The radial evolution of $\sigma _{c}$ and $\sigma _{r}$ on the circle plot of figure 1 provides simple, persuasive evidence that our reflection-turbulence model captures key aspects of solar-wind evolution. Observations robustly show that solar-wind fluctuations are concentrated near the circle's bottom-right quadrant edge, evolving from $(\sigma _c,\sigma _r) =(1,0)$ close to the sun towards $(\sigma _c,\sigma _r) =(0,-1)$ at large radii (Bavassano, Pietropaolo & Bruno Reference Bavassano, Pietropaolo and Bruno1998; Bruno et al. Reference Bruno, D'Amicis, Bavassano, Carbone and Sorriso-Valvo2007; Bruno & Carbone Reference Bruno and Carbone2013; Wicks et al. Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013). This behaviour agrees with our simulations and phenomenological arguments (figure 1): during its imbalanced phase, the system remains close to the circle's edge because reflection generates $\boldsymbol {z}^{-}$ fluctuations that are anti-aligned with their $\boldsymbol {z}^+$ source (negative $\sigma _{\theta }$; see (3.3)), evolving into the $\sigma _{r}\simeq -1$ state at late times as the long-wavelength expansion-dominated modes start dominating the dynamics (§ 5). In addition, faster wind is observed to be concentrated near the middle right (large $\sigma _{c}$, small $\sigma _{r}$), while slower streams are concentrated near the bottom (large $\sigma _{r}$, small $\sigma _{c}$; Bavassano et al. Reference Bavassano, Pietropaolo and Bruno1998; Bruno et al. Reference Bruno, D'Amicis, Bavassano, Carbone and Sorriso-Valvo2007), which fits straightforwardly into the idea described above that fast-wind streams have larger $\chi _{\rm exp}$, thus spending longer in the imbalanced phase. While we are, again, not the first to speculate on the relevance of expansion to these observations (Zhou & Matthaeus Reference Zhou and Matthaeus1990; Oughton & Matthaeus Reference Oughton and Matthaeus1995; Hollweg & Isenberg Reference Hollweg and Isenberg2007; Grappin et al. Reference Grappin, Verdini and Müller2022), we believe ours are the first simulations to highlight this feature, particularly the evolution into the balanced $\sigma _{r}\simeq -1$ state and the importance of the Elsässer alignment $\sigma _\theta$. Note that one can take advantage of the high alignment observed to deduce a radial evolution law for the normalized residual energy. Using the fact that the system remains close to the circle's edge where $\sigma _r^2 \sim 1-\sigma _c^2$, one finds that

(7.3)\begin{equation} \sigma_{r} \sim \begin{cases} -\dfrac{2\chi_{\rm exp}^{{-}2}}{1+\chi_{\rm exp}^{{-}2}}, & \chi_{\rm exp}\gtrsim 1, \\[10pt] -1 , & \chi_{\rm exp} \lesssim 1. \end{cases} \end{equation}

As for the relation (7.2), this prediction should be directly testable with spacecraft observations.

7.2.5 Solar-wind heating

Our study also has application to the understanding of solar-wind heating, although more work is needed. In fast-wind streams the observed radial decrease of the proton temperature $T$ is slower than in adiabatic cooling, indicating that the plasma is heated as it moves outward from the Sun (Marsch et al. Reference Marsch, Mühlhäuser, Schwenn, Rosenbauer, Pilipp and Neubauer1982; Freeman Reference Freeman1988; Tu Reference Tu1988), presumably by the dissipation of fluctuations (Verma, Roberts & Goldstein Reference Verma, Roberts and Goldstein1995; Vasquez et al. Reference Vasquez, Smith, Hamilton, MacBride and Leamon2007; Chen et al. Reference Chen2020; Halekas et al. Reference Halekas2023). Such turbulent heating appears to be less important in slower-wind streams, although results remain controversial (Freeman Reference Freeman1988; Totten, Freeman & Arya Reference Totten, Freeman and Arya1995; Hellinger et al. Reference Hellinger, Matteini, Štverák, Trávníček and Marsch2011). Within the RMHD EBM model, the heating rate per unit volume is $Q = -\rho \langle \boldsymbol {z}^+\boldsymbol {\cdot }\boldsymbol {D}^++\boldsymbol {z}^{-} \boldsymbol {\cdot }\boldsymbol {D}^{-}\rangle /2$, where $\boldsymbol {D}^{\pm }$ represents the hyperviscous terms included to dissipate energy at small scales (2.10). Converting to wave-action variables and using the total energy conservation (2.8) (see also Perez et al. Reference Perez, Chandran, Klein and Martinović2021) gives

(7.4)\begin{equation} Q ={-} \rho\frac{\dot{a}}{a}\left( \frac{\partial \tilde{E}}{\partial a} + \frac{\tilde{E}^{r}}{a}\right).\end{equation}

During the imbalanced phase at high $\chi _{\rm exp}$, when $\tilde {E}^{r}\ll \tilde {E}\approx \tilde {E}^+$, the phenomenology of § 4.1 predicts $\tilde {E}^+\propto a^{-1}$ and $\rho \propto a^{-2}$, so that $Q\propto a^{-5}$. Then, as the system transitions into the balanced phase, the heating rate drops significantly as the system becomes dominated by slowly evolving Alfvén vortices. At late times we measure a small residual nonlinear dissipation that causes $\tilde {E}\propto a^{0.8}$ (rather than the $\tilde {E}\propto a^{1}$ predicted by linear theory), implying a heating rate that flattens to $Q\propto a^{-3.2}$.

Whether these results are consistent with observations remains unclear. Most observational studies have inferred heating rates by fitting power-law profiles to the observed temperatures, then comparing the inferred scalings to ‘adiabatic’ profiles that would occur in the absence of heating: $T\propto R^{-4/3}$ for an isotropic fluid (i.e. if the perpendicular and parallel temperatures are well coupled, $T_\perp \sim T_\|$), or $T_\perp \propto R^{-2}$ for a collisionless plasma (or more generally, $T_\perp \propto B$). A $Q\propto a^{-\alpha }$ heating profile with $\alpha \geqslant 5$ ($\alpha \geqslant 13/3$ for an isotropic fluid) is too steep to lead to a power-law temperature profile that differs from the adiabatic profile at asymptotically large $R$; however, depending on the magnitude of $Q$, almost any local scaling of $T$ can be realized (it just does not vary as a power law over a wide range in $R$). This, combined with the effects of averaging over different streams with different $\chi _{\rm exp}$, makes it unclear whether the difference between the high-$\chi _{\rm exp}$ prediction ($Q\propto a^{-5}$) and the classic result that $\alpha \approx 3.8\rightarrow 4$ (Freeman Reference Freeman1988; Totten et al. Reference Totten, Freeman and Arya1995; Hellinger et al. Reference Hellinger, Matteini, Štverák, Trávníček and Marsch2011) should signal the importance of other physics, or not. Indeed, more complex models based on a similar phenomenology (Cranmer Reference Cranmer2009; Chandran et al. Reference Chandran, Dennis, Quataert and Bale2011) reproduce observed temperature profiles reasonably well out to ${\simeq }1\ {\rm AU}$, and Montagud-Camps et al. (Reference Montagud-Camps, Grappin and Verdini2020) found that expanding-box turbulence with $\chi _{\rm exp0}\simeq 5$ and a high enough Mach number can reasonably reproduce the observed $T\propto R^{-1}$ temperature profile over a factor of several in radius. Also of interest is the transition around $\chi _{\rm exp}\sim 1$, where we predict that the decrease in heating rate with radius should slow to eventually approach $Q\propto a^{-3}$ as the heating stops. If slow-wind streams have smaller $\chi _{\rm exp}$ as suggested above, the general trend could be consistent with the observation of closer-to-adiabatic evolution in slow wind (a flatter power-law profile of $Q$ will not be measurable if its magnitude is too small), as well as recent measurements showing the much greater importance of wave heating in fast, compared with slow streams (Halekas et al. Reference Halekas2023).

Overall, while plausibly consistent, more work is needed, particularly to quantify the relevance reflection-driven turbulence compared with other effects. Those of particular relevance for further study include pickup ions, which are thought to contribute significantly at larger radii (Gazis et al. Reference Gazis, Barnes, Mihalov and Lazarus1994; Zank et al. Reference Zank, Adhikari, Zhao, Mostafavi, Zirnstein and McComas2018), stream interactions across a range of radii (e.g. Roberts et al. Reference Roberts, Goldstein, Matthaeus and Ghosh1992; Breech et al. Reference Breech, Matthaeus, Minnie, Bieber, Oughton, Smith and Isenberg2008) and parametric decay in highly imbalanced regions (Chandran et al. Reference Chandran, Foucart and Tchekhovskoy2018; Shoda et al. Reference Shoda, Suzuki, Asgari-Targhi and Yokoyama2019).