Reflection-driven turbulence in the super-Alfvénic solar wind

Abstract

In magnetized, stratified environments such as the Sun's corona and solar wind, Alfvénic fluctuations ‘reflect’ from background gradients, enabling nonlinear interactions that allow their energy to dissipate into heat. This process, termed ‘reflection-driven turbulence’, likely plays a key role in coronal heating and solar-wind acceleration, explaining a range of detailed observational correlations and constraints. Building on previous works focused on the inner heliosphere, here we study the basic physics of reflection-driven turbulence using reduced magnetohydrodynamics in an expanding box – the simplest model that can capture local turbulent plasma dynamics in the super-Alfvénic solar wind. Although idealized, our high-resolution simulations and simple theory reveal a rich phenomenology that is consistent with a diverse range of observations. Outwards-propagating fluctuations, which initially have high imbalance (high cross-helicity), decay nonlinearly to heat the plasma, becoming more balanced and magnetically dominated. Despite the high imbalance, the turbulence is strong because Elsässer collisions are suppressed by reflection, leading to ‘anomalous coherence’ between the two Elsässer fields. This coherence, together with linear effects, causes the growth of ‘anastrophy’ (squared magnetic potential) as the turbulence decays, forcing the energy to rush to larger scales and forming a ‘$1/f$-range’ energy spectrum in the process. Eventually, expansion overcomes the nonlinear and Alfvénic physics, forming isolated, magnetically dominated ‘Alfvén vortices’ with minimal nonlinear dissipation. These results can plausibly explain the observed radial and wind-speed dependence of turbulence imbalance (cross-helicity), residual energy, fluctuation amplitudes, plasma heating and fluctuation spectra, as well as making a variety of testable predictions for future observations.

1 Introduction

The mechanisms that heat and accelerate the solar wind remain mysterious, or at least controversial (Cranmer & Winebarger 2019). In order to explain decades of in-situ spacecraft data, particularly local temperature measurements and the high speeds of fast-wind streams, there must exist an energy source to heat the plasma even at large distances from the solar surface. A leading paradigm for explaining this extended heating is Alfvénic turbulence, in which the energy is provided by Alfvén waves launched from the low solar atmosphere by photospheric motions or magnetic reconnection (Axford & McKenzie 1992; De Pontieu et al. 2007). As these waves propagate outwards, away from the Sun, they become turbulent, causing their energy to cascade to smaller scales and dissipate (Velli, Grappin & Mangeney 1989; Dmitruk et al. 2002; Dmitruk & Matthaeus 2003; Cranmer & van Ballegooijen 2005). The resulting turbulent heating increases the plasma pressure, which, along with the wave pressure, accelerates the solar wind away from the Sun (Tu 1987, 1988; Cranmer, van Ballegooijen & Edgar 2007; Verdini et al. 2010).

Although plausible, particularly given the extended turbulent-like fluctuations observed in the solar-wind plasma (Belcher & Davis 1971; Bruno & Carbone 2013; Kiyani, Osman & Chapman 2015; Chen 2016), a particular difficulty with this model lies in the robustness of Alfvénic fluctuations: in a homogenous plasma, Alfvén waves propagating in the same direction do not interact with one another or damp out, even at large amplitudes and/or when their wavelength is well below the mean free path (Barnes & Hollweg 1974; Kulsrud 1983). Turbulence, as likely needed to dissipate their energy, thus arises only via interactions between the two ‘Elsässer’ fields $\boldsymbol {z}^\pm$, which are the counter-propagating linear eigenmodes in a homogenous plasma (Elsasser 1950; Iroshnikov 1963; Kraichnan 1965). With the Sun supplying energy only in outwards-propagating waves dominated by one Elsässer field, some source of the other Elsässer field is needed to generate turbulence that could explain the observed heating. One possible mechanism for enabling this process is reflection arising from the radial variation in the background Alfvén speed $v_{A}$ (Heinemann & Olbert 1980; Zhou & Matthaeus 1989; Velli 1993). The turbulence that results due to this interaction between outwards and reflected waves is generally referred to as ‘reflection-driven turbulence’ (Velli et al. 1989). Phenomenological models and simulations suggest that the paradigm can broadly explain many observed local and global features of the solar wind (e.g. Dmitruk et al. 2002; Suzuki & Inutsuka 2005; Cranmer & van Ballegooijen 2005; Verdini & Velli 2007; Chandran & Hollweg 2009; Verdini et al. 2012; Perez & Chandran 2013; van der Holst et al. 2014; van Ballegooijen & Asgari-Targhi 2017), although there remain important unresolved issues and other physical effects of importance (e.g. Lionello et al. 2014; Zank et al. 2017; Shoda et al. 2019; Asgari-Targhi et al. 2021; Chandran 2021). Similar mechanisms may also play a key role in other astrophysical systems with large density gradients and strong magnetic fields, particularly compact-object accretion flows, which are known to possess hot, compact corona that are likely fed by strong fluctuations in the disk below (Reis & Miller 2013; Chandran, Foucart & Tchekhovskoy 2018).

The goal of this work is to study reflection-driven turbulence from the simplest standpoint possible, elucidating its key features in a minimally complex setting. This differs from previous studies, which have usually used either phenomenological turbulence/transport models (e.g. Zank, Matthaeus & Smith 1996; Cranmer & van Ballegooijen 2005; Suzuki & Inutsuka 2005; Verdini & Velli 2007; van der Holst et al. 2014; Réville et al. 2020) or radially extended, inhomogenous ‘flux-tube’ simulations (e.g. Dmitruk et al. 2002; Dmitruk, Gómez & Matthaeus 2003; Perez & Chandran 2013; van Ballegooijen & Asgari-Targhi 2016) to attempt to realistically match observed parameters and regimes of the corona and solar wind. Both perspectives – the simplified and the realistic – are useful, but the former has seen less attention in previous literature. Focusing on the simplified perspective is especially relevant because reflection-driven turbulence is neither purely decaying nor forced (the two limits usually considered in turbulence studies), but as we shall see, involves features of both limits at the same time. This means that care is needed when applying intuitions and ideas from broader turbulence research. While many different effects are undeniably important in a system as complex as the solar wind, we argue the process of neglecting physics – even that which might be important – can be crucial for uncovering interesting effects that could be missed in more complete models. To put this work in context, the simple phenomenological picture put forward to explain our simulations contains only Alfvénic and non-Wentzel–Kramers–Brillouin (WKB) reflection physics, similar to previous flux-tube reflection models (e.g. Suzuki & Inutsuka 2005; Verdini & Velli 2007) but with modified predictions at large radius because the system approaches nonlinear solutions and halts the turbulent cascade. Another collection of literature (e.g. Zank et al. 1996; Matthaeus et al. 1999; Breech et al. 2008) considers the effects of ‘mixing, expansion, compression and shear’ (MECS) on turbulent transport; our work here contains only the mixing and expansion from the spherical flow, neglecting the effects of large-scale stream–stream shear or compression. While such effects may sometimes be of importance, we will argue they are not needed to explain many basic features, including the evolution of the cross-helicity, residual energy and fluctuation amplitude.

Our numerical approach is to use the so-called ‘expanding-box model’ (EBM) (Grappin, Velli & Mangeney 1993), which tracks a small parcel of plasma as it flows away from the Sun. The version of the EBM we use applies to regions beyond the Alfvén radius (or surface) $R_{A}$ where the solar-wind speed $U$ overtakes the Alfvén speed and becomes approximately constant with radius $R$. This local EBM approach differs from most previous work on reflection-driven turbulence although the effect of expansion on turbulence has been studied with the EBM in various contexts; (e.g. Grappin & Velli 1996; Dong, Verdini & Grappin 2014; Montagud-Camps, Grappin & Verdini 2018; Grappin, Verdini & Müller 2022; Johnston et al. 2022). A disadvantage of the EBM is that our results cannot be applied directly to the solar-wind acceleration region (although some aspects may prove translatable); an advantage is the simplicity of using a homogenous, periodic domain, which allows for much higher numerical resolutions and decreases the number of free parameters while capturing many of the essential physical ingredients. In addition, our results seem to explain a variety of disparate observations from in-situ spacecraft measurements at $R>R_{A}$, some of which have been missed in previous theoretical works because of the focus on lower-altitude acceleration regions. We argue that these observational comparisons provide persuasive evidence that reflection-driven turbulence controls important aspects of solar-wind turbulent evolution beyond $R_{A}$, as well as providing a number of testable and falsifiable predictions for future works.

As well as the contributions described above, our main novel result is that reflection-driven turbulence precipitates a strong inverse energy transfer as it decays. This feature, which we argue is a consequence of an anomalous¹ conservation law for the squared parallel magnetic vector potential (‘anastrophy’), causes initially small-scale outwards-propagating fluctuations to rush to large scales as they decay, forming a ${\propto }k_{\perp }^{-1}$ spectrum in the process (here $k_{\perp }$ is the wavenumber perpendicular to the background magnetic field). This suggests the observed large-scale fluctuations that dominate the solar-wind turbulence spectrum can develop in situ as the wind propagates, which may be important if low-frequency waves are unable to effectively propagate through the chromosphere-coronal transition due to large local gradients in the Alfvén speed (Leroy 1981; Velli 1993; van Ballegooijen & Asgari-Targhi 2017; Réville, Tenerani & Velli 2018). Another new result concerns the asymptotic evolution of the turbulence at large radii, where it becomes governed by large-scale magnetically dominated ‘Alfvén vortices’ (Petviashvili & Pokhotelov 1992; Alexandrova 2008). These structures, which are approximate nonlinear solutions and so dissipate into heat only very slowly, tend to freeze into the plasma at late times, growing continuously as the plasma expands in a way that is consistent with the observed slow decay of turbulent fluctuation amplitudes at large radii (e.g. Zank et al. 1996).

The remainder of the paper is organized as follows. Section 2 describes the basic expanding-box reduced magnetohydrodynamic (RMHD) model that we use throughout this work. We outline the useful ‘wave-action’ form (§ 2.1.1), which facilitates analysis by factoring out the linear WKB wave evolution brought about by expansion (though non-WKB physics is still retained in the model). We then explain the numerical method, key parameters of the system and the initial conditions used for the simulations. Section 3 then presents a brief overview of how the turbulence evolves, focusing on globally averaged quantities such as the energy, imbalance (normalized cross-helicity) and residual energy. We will see that the evolution splits into two distinct phases, evolving from one nonlinear solution of homogenous MHD (pure outwards-propagating waves, high imbalance) to another (magnetically dominated Alfvén vortices). In § 4 we examine the imbalanced phase, starting with a simple phenomenology based on previous works (Dmitruk et al. 2002; Verdini & Velli 2007; Chandran & Hollweg 2009) to understand the observed dynamics. We compare these phenomenological ideas to the simulations’ time evolution (§ 4.1), spectra (§ 4.2) and frequency spectra (§ 4.3), diagnosing how the suppression of wave collisions leads to ‘anomalous coherence’, enabling strong turbulence despite the high imbalance. In § 4.4 we then examine the inverse transfer in detail, presenting a theoretical argument based on anastrophy to explain the observed results. The balanced, magnetically dominated phase is examined in § 5, starting with a focus on linear expansion-dominated (long-wavelength, non-WKB) physics (§ 5.1). This linear physics controls the late-stage evolution of the system because the system self-organizes to minimize its nonlinearity, explaining the strong dominance of magnetic over kinetic energy and various other features of its evolution (as well as a number of solar-wind observations). That this system does indeed morph into nonlinear solutions is proved numerically (and argued theoretically) by directly fitting structures that grow in the simulation (§ 5.3).

The paper contains a lot of detail about various aspects of the evolution. Therefore, in § 6 we provide a summary of how this work relates to previous literature on solar-wind turbulence, followed by an extended discussion of the observational relevance of our findings in § 7. The latter covers explanations of various existing observational results, such as the observed radial evolution and wind-speed dependence of imbalance and residual energy, as well as making predictions that can be tested in future works to better understand the successes and limitations of the reflection-driven turbulence model. We conclude in § 8.

2 Methods

2.1 The expanding RMHD model

We wish to describe the turbulent dynamics of a plasma advected by an expanding wind and threaded by a mean magnetic field $\bar {\boldsymbol {B}}$ using the simplest possible formalism. We therefore assume that $\bar {\boldsymbol {B}}$ is radial, and that the fluctuations in the total field $\boldsymbol {B}$ and plasma velocity $\boldsymbol {u}$ are transverse and non-compressive, with characteristic scales well above the ion gyroscale (i.e. the fluctuations are polarized like shear-Alfvén waves). We assume that the mean expanding flow of the wind $\boldsymbol {U}$ is also radial, constant and much larger than the Alfvén speed $v_{A}\equiv |\bar {\boldsymbol {B}}|/\sqrt {4{\rm \pi} \rho }$, where $\rho$ is the mass density of the plasma. These assumptions about $\boldsymbol {u}$ and $\boldsymbol {B}$ apply reasonably well to the solar-wind plasma in regions with $\mathcal {M}_{A} \equiv |\boldsymbol {U}|/v_{A}\gtrsim 1$ (i.e. beyond the Alfvén point) and where the Parker spiral is still well aligned with the radial direction (Parker 1965). Even with such simplifications, simulating such dynamics using an absolute frame of reference and over a large radial distance remains extremely costly in terms of computer power (Perez & Chandran 2013; van Ballegooijen & Asgari-Targhi 2016; Chandran & Perez 2019). We circumvent this difficulty by considering the turbulence dynamics in a frame co-moving with the spherically expanding flow – the so-called EBM (Grappin et al. 1993). Assuming that the domain is small compared with the heliocentric distance, the curvature of surfaces perpendicular to the radially expanding flow can be neglected, allowing the use of Cartesian coordinates and periodic boundary conditions in all three directions. The resulting savings in numerical cost are redeployed to resolve the turbulence across a range of scales of unprecedented breadth.

These approximations lead to equations that take a form similar to RMHD (Kadomtsev & Pogutse 1973; Schekochihin et al. 2009), with two modifications. First, there appear additional linear terms proportional to $\boldsymbol {U}_{\perp }$, which is the compressive part of the mean radial velocity perpendicular to the radial direction at the centreline of the simulation domain, which acts to expand the domain as it moves outwards (note that the non-radial part of $\bar {\boldsymbol {B}}$ can be neglected because $|\boldsymbol {U}|\gg v_{A}$ and due to the small spatial domain). Second, the perpendicular gradient operator is modified to account for the increasing lateral stretching of the plasma with distance: $\hat {\boldsymbol {\nabla }} \equiv (a^{-1}\partial _x,a^{-1}\partial _y, \partial _z)$, where we use the local-box spatial coordinates $(x,y,z)$ and align the $z$ axis with the outwards radial direction at the centreline of the simulation domain. Here $a$ is defined as the heliospheric distance $R$ of the co-moving frame, normalized by the initial radial distance $R_0$ (equivalently, it is the perpendicular size of the domain):

(2.1)\begin{equation} a(t)=\dfrac{R(t)}{R_{0}}=\dfrac{R_{0}+Ut}{R_{0}}=1+\dot{a}t.\end{equation}

Here $\dot {a}=\partial a/\partial t = U/R_0$ is a constant for constant $\boldsymbol {U}$. Noting that $\boldsymbol {U}_\perp = (\dot {a}/a)(x\hat {\boldsymbol {x}}+y\hat {\boldsymbol {y}})$, one finds that the magnetic field, $\boldsymbol {B} =\bar {\boldsymbol {B}}+\boldsymbol {B}_{\perp }=B_z\hat {\boldsymbol {z}}+\boldsymbol {B}_{\perp }$, and the part of the perpendicular flow velocity that remains after the Galilean transformation, $\boldsymbol {u}_\perp = \boldsymbol {u} - \boldsymbol {U}$, evolve as (Grappin et al. 1993)

(2.2)\begin{align} \dfrac{{\rm d}\boldsymbol{u}_{{\perp}}}{{\rm d} t} + \dfrac{\hat{\boldsymbol{\nabla}}_{{\perp}} p}{\rho}-\dfrac{\boldsymbol{B}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}\boldsymbol{B}_\perp}{4{\rm \pi} \rho}& ={-}\boldsymbol{u}_{{\perp}}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}_{{\perp}} \boldsymbol{U}_{{\perp}}\nonumber\\ & ={-}\dfrac{\dot{a}}{a}\boldsymbol{u}_\perp, \end{align}

(2.3)\begin{align} \dfrac{{\rm d}\boldsymbol{B}}{{\rm d} t}-\boldsymbol{B}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}_{{\perp}}\boldsymbol{u}_{{\perp}} & ={-}\boldsymbol{B}\hat{\boldsymbol{\nabla}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{U}_{{\perp}} +\boldsymbol{B}_{{\perp}}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}_{{\perp}} \boldsymbol{U}_{{\perp}} \nonumber\\ & ={-} 2\dfrac{\dot{a}}{a} B_z \hat{\boldsymbol{z}} -\dfrac{\dot{a}}{a} \boldsymbol{B}_\perp, \end{align}

where ${\rm d}/{\rm d} t = \partial /\partial t + \boldsymbol {u}_{\perp }\boldsymbol {\cdot }\hat {\boldsymbol {\nabla }}_{\perp }$. The total pressure $p$, which includes both magnetic and thermal pressures, cancels the compressive part of the nonlinear terms to enforce the incompressibility of the motions $\hat {\boldsymbol {\nabla }}_{\perp } \boldsymbol {\cdot }\boldsymbol {u}_\perp =0$ (Schekochihin et al. 2009). Defining the subscript $0$ to refer to a quantity at $t=0$ ($a=1$), conservation of mass and magnetic flux imply that $\rho =\rho _0/a^2$ and $B_z=B_{z0}/a^2$ (the latter being the solution of the $\hat {\boldsymbol {z}}$ component of (2.3)), so that $v_{A} =v_{A0}/a$ (Grappin et al. 1993). Note that because $\rho =\rho _0/a^2$, the perpendicular friction-like term in (2.3) associated with the spherical expansion ($-\dot {a}/a \boldsymbol {B}_\perp$) vanishes if one instead expresses the perpendicular magnetic field in velocity units $\boldsymbol {b}_{\perp }=\boldsymbol {B}_{\perp }/\sqrt {4{\rm \pi} \rho }$ using $\partial _t \boldsymbol {b}_{\perp }=(4{\rm \pi} \rho )^{-1/2}\partial _t \boldsymbol {B}_{\perp }+\dot {a}/a\, \boldsymbol {b}_{\perp }$. Because $\boldsymbol {u}_\perp$ is damped via $-\dot {a}/a\, \boldsymbol {u}_\perp$, this produces differential damping of the perpendicular magnetic and kinetic fluctuations during the radial transport.

The most important impact of expansion is that it causes Alfvénic reflection. This can be seen by considering the Elsässer variables $\boldsymbol {z}^{\pm }=\boldsymbol {u}_{\perp }\pm \boldsymbol {b}_{\perp }$, which evolve as

(2.4)\begin{equation} \dfrac{\partial \boldsymbol{z}^{{\pm}}_{{\perp}}}{\partial t} \pm v_{A}\dfrac{\partial\boldsymbol{z}^{{\pm}}_{{\perp}}}{\partial z}+ \boldsymbol{z}^{{\mp}}_{{\perp}} \boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}_{{\perp}}\boldsymbol{z}^{{\pm}}_{{\perp}}+\dfrac{\hat{\boldsymbol{\nabla}}_{{\perp}} p}{\rho}={-}\dfrac{1}{2}\dfrac{\dot{a}}{a}\left( \boldsymbol{z}^+_{{\perp}}+\boldsymbol{z}^{-}_{{\perp}}\right). \end{equation}

We have taken $\bar {\boldsymbol {B}}$ to point in the negative radial direction ($B_z<0$ with $v_{A}= |B_z|/\sqrt {4{\rm \pi} \rho }$), so that $\boldsymbol {z}^+_{\perp }$ perturbations propagate outwards in the absence of reflection.² We see that the additional linear terms proportional to $\boldsymbol {U}_{\perp }$ appearing in (2.3) and (2.2) couple $\boldsymbol {z}^+_{\perp }$ and $\boldsymbol {z}^{-}_{\perp }$ perturbations through the final term in (2.4), with important consequences for their nonlinear evolution.

The RMHD-like form of (2.4) can be derived from compressible MHD, or even the full Vlasov–Maxwell system, by considering anisotropic ($\partial _z\ll \boldsymbol {\nabla }_\perp$) and small-amplitude fluctuations, even when the compressive fluctuations have similar amplitudes to the Alfvénic ones (Schekochihin et al. 2009; Oughton, Matthaeus & Dmitruk 2017). However, observations of large-amplitude field reversals, or switchbacks, show that highly imbalanced turbulence in the near-Sun solar wind is often large amplitude ($\delta B/B\sim 1$), violating the latter assumption, although such structures do seem to be anisotropic (Horbury et al. 2020; indeed, a spherically polarized fluctuation must be anisotropic to allow a field reversal $\delta B_\| \sim \bar {B}$; Mallet et al. 2021). The reduced nature of the model is, therefore, a priori unjustified for describing the super-Alfvénic wind where fluctuation amplitudes are large and the Parker spiral becomes significant. Our motivation to use this model despite these deficiencies stems from the hope that, by isolating the Alfvénic dynamics from other physics, RMHD can provide an approximate description of key dynamical effects while remaining conceptually simple (see Dmitruk, Matthaeus & Oughton (2005) for a comparison of compressible and RMHD in simulations). While the comparison to some previous EBM results (§ 6) and observations (§ 7) provides an indication that this hope could be warranted, clearly more work with a more complete model is needed. In § 8.2 we examine how the simplifying assumptions underlying this model may influence the conclusions drawn in this paper.

2.1.1 Wave-action form

It is convenient to rewrite (2.4) in terms of the so-called wave-action Elsässer variables (Heinemann & Olbert 1980), defined as

(2.5)\begin{equation} \tilde{\boldsymbol{z}}^{{\pm}} \doteq a^{1/2}\boldsymbol{z}^{{\pm}}_\perp{\propto} \frac{ \boldsymbol{z}^{{\pm}}_\perp }{\sqrt{\omega_{A}}} , \end{equation}

where $\omega _{A}=k_{z}v_{A}$ is the Alfvén frequency of a mode of wavenumber $k_{z}$. The second expression emphasizes the relationship to the wave-action density (Whitham 1965), which is $|\boldsymbol {z}^{\pm }|^{2}/\omega _{A}$ for a population of $\boldsymbol {z}^{\pm }$ fluctuations at some $k_{z}$, and is conserved in the limit of high-frequency/short-wavelength waves (but note that we do not neglect the long-wavelength non-WKB physics in the change of variables). This highlights how the extra $a^{-1/2}$ factor compensates the decay of the $\boldsymbol {z}^{\pm }_{\perp }$ that arises because of the decreasing Alfvén frequency as the system expands, making $\tilde {\boldsymbol {z}}^{\pm }$ the natural variables in which to consider turbulent-decay dynamics. Equation (2.4) then take the form

(2.6)\begin{equation} \dot{a}\dfrac{\partial \tilde{\boldsymbol{z}}^{{\pm}}}{\partial a}\pm v_{A} \dfrac{\partial \tilde{\boldsymbol{z}}^{{\pm}}}{\partial z}+\dfrac{1}{a^{1/2}}\left(\tilde{\boldsymbol{z}}^{{\mp}}\boldsymbol{\cdot}\hat{\boldsymbol{\nabla}}_{{\perp}} \tilde{\boldsymbol{z}}^{{\pm}} + \dfrac{\hat{\boldsymbol{\nabla}}_{{\perp}} p}{\rho} \right)={-}\dfrac{\dot{a}}{2a}\tilde{\boldsymbol{z}}^{{\mp}}. \end{equation}

These equations can be equivalently derived from the ‘flux-tube’ RMHD equations used by Perez & Chandran (2013), Chandran & Perez (2019) (see also Verdini & Velli 2007; Chandran & Hollweg 2009) by identifying their $\boldsymbol {g}$ and $\boldsymbol {f}$ with $\tilde {\boldsymbol {z}}^+$ and $\tilde {\boldsymbol {z}}^{-}$, respectively, assuming $v_{A}\ll U$ and $\eta \doteq \rho /\rho |_{U=v_{A}}\ll 1$, and converting $\dot {a}\partial /\partial a$ in (2.6) into the derivative in the stationary frame $\partial /\partial t + U \partial /\partial R$.

For the remainder of this paper we will usually use wave-action variables with lengths and gradients defined in the co-moving frame, which does not change with $a$. With this in mind, it is sometimes helpful to explicitly expand $\hat {\boldsymbol {\nabla }}_{\perp }=a^{-1}\tilde {\boldsymbol {\nabla }}_\perp$ and $v_{A}=v_{A0}/a$, in order to remove the hidden $a$ dependence of these terms in (2.6). Written in terms of $\ln a$, (2.6) takes a form that is similar to standard RMHD in a fixed-size domain, but with reflection terms and a time-variable coefficient $a^{-1/2} = {\rm e}^{-\ln a/2}$ multiplying the nonlinear term

(2.7)\begin{equation} \dot{a}\dfrac{\partial \tilde{\boldsymbol{z}}^{{\pm}}}{\partial \ln a}\pm {v_{A0}} \dfrac{\partial \tilde{\boldsymbol{z}}^{{\pm}}}{\partial z}+\frac{1}{a^{1/2}}\tilde{\boldsymbol{z}}^{{\mp}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}_\perp \tilde{\boldsymbol{z}}^{{\pm}} + \tilde{\boldsymbol{\nabla}}_\perp \tilde{p}={-}\frac{\dot{a}}{2}\tilde{\boldsymbol{z}}^{{\mp}}, \end{equation}

where $\tilde {p}=p/\rho$ enforces $\tilde {\boldsymbol {\nabla }}_\perp \boldsymbol {\cdot }\tilde {\boldsymbol {z}}^{\mp }=0$. It is often helpful to consider the turbulent evolution from the perspective of (2.7), multiplying lengths by $a$ and using (2.5) to convert back to physical quantities as need be. We similarly define wave-action velocities and magnetic fields, $\tilde {\boldsymbol {u}}_{\perp } = a^{1/2} \boldsymbol {u}_{\perp }$ and $\tilde {\boldsymbol {b}}_{\perp } = a^{1/2}\boldsymbol {b}_{\perp }$, respectively.

Throughout this paper we use the tilde $\tilde {\cdot }$ to denote both wave-action-normalized fields and length scales defined in the co-moving frame (like $\tilde {\boldsymbol {\nabla }}_\perp$). Because we have not transformed time in deriving (2.6) or (2.7), time scales and frequencies are not adorned with a tilde, and can be equivalently defined in either the co-moving or physical frame with either wave-action or physical variables, as convenient. The same is true for dimensionless quantities and parallel length scales.

2.1.2 Conserved quantities

Unlike homogeneous RMHD, individual wave-action Elsässer energies $\tilde {E}^{\pm } \equiv \langle |\tilde {\boldsymbol {z}}^{\pm } |^2 \rangle /4$ are not conserved in the presence of expansion. (Here and in the following, angle brackets $\langle \dots \rangle$ denote a volume average over the expanding box in the co-moving frame.) The reflection terms can act as a source or a sink of wave-action energy, depending on the sign of the correlation between the Elsässer fields or residual energy $\tilde {E}^r =\langle \tilde {\boldsymbol {z}}^+\boldsymbol {\cdot }\tilde {\boldsymbol {z}}^{-}\rangle /2 = \tilde {E}^{u} - \tilde {E}^{b}$ (we define also the wave-action kinetic and magnetic energies, $\tilde {E}^{u} = \langle |\tilde {\boldsymbol {u}}_{\perp }|^{2}\rangle /2$ and $\tilde {E}^{b} = \langle |\tilde {\boldsymbol {b}}_{\perp }|^{2}\rangle /2$, respectively). Specifically, one finds from (2.6),

(2.8)\begin{equation} \dot{a}\dfrac{\partial \tilde{E}^{{\pm}}}{\partial a} ={-}\dfrac{\dot{a}}{4a}\langle \tilde{\boldsymbol{z}}^+\boldsymbol{\cdot}\tilde{\boldsymbol{z}}^{-}\rangle={-}\dfrac{\dot{a}}{2a}\tilde{E}^r.\end{equation}

For the same reason, the total energy $\tilde {E} = \tilde {E}^++\tilde {E}^{-}=\tilde {E}^{u}+\tilde {E}^{b}$ is not conserved. In contrast, one sees that the reflection sources cancel out for the wave-action cross-helicity $\tilde {E}^{c}= \tilde {E}^+-\tilde {E}^{-} = \langle \tilde {\boldsymbol {u}}_{\perp }\boldsymbol {\cdot }\tilde {\boldsymbol {b}}_{\perp }\rangle$, which therefore remains, as in the homogeneous case, an ideal invariant (Verdini & Velli 2008),

(2.9)\begin{equation} \dfrac{\partial \tilde{E}^{c}}{\partial a} =0.\end{equation}

We note that although the fluctuation energy is not conserved, one can show using the full system of equations (without making the expanding-box approximation) that the energy gained or lost by the fluctuations is compensated by an equal and opposite change in the energy of the background flow (Chandran, Schekochihin & Mallet 2015; Perez et al. 2021).

2.2 Numerical method and set-up

Taking advantage of the periodic boundary conditions, we solve (2.2) and (2.3) (or equivalently, (2.4), (2.6) or (2.7)) with a modified version of the Fourier pseudo-spectral code TURBO (Teaca et al. 2009). We solve (2.2) and (2.3) in potential form, using time-dependent $k_\perp$ to account for the expansion and advancing in time with a third-order modified Williamson algorithm (a four-step, low-storage Runge–Kutta method; Williamson 1980) for the nonlinear terms and implicitly evaluate the linear terms exactly. The simulation domain is a cube of size $L_\perp =L_z=2{\rm \pi}$ with a resolution $n_{\perp }^{2}\times n_{z}$. Note that the system (2.6) has a rescaling symmetry, whereby all relative fluctuation amplitudes can be arbitrarily rescaled as long as the ratios of all perpendicular to parallel scales are rescaled by the same amount. Therefore, the parallel and perpendicular units of length are independent. The code units are set by this and by $v_{A0} = 2{\rm \pi}$. Nonlinear terms are partially dealiased using a phase-shift method (Patterson & Orszag 1971). The main simulations presented below will use a spatial resolution of $n_{\perp }^{2}\times n_{z}=1536^{2}\times 128$ for the full simulation evolution, but are refined to $n_{\perp }^{2}\times n_{z}=8192^{2}\times 256$ around specified radii (time) of interest and allowed to evolve briefly, in order to resolve spectra at smaller scales.

We add a form of dissipation (‘hyperviscosity’),

(2.10)\begin{equation} -\nu^{{\pm}}_{{\perp}}\hat{\boldsymbol{\nabla}}_{{\perp}}^{6}\tilde{\boldsymbol{z}}^{{\pm}} -\nu^{{\pm}}_{z}\partial_{z}^{6}\tilde{\boldsymbol{z}}^{{\pm}}, \end{equation}

to the right-hand side of (2.6) to absorb the turbulent energy at small scales. The hyperviscosity coefficients $\nu ^{\pm }_{\perp }$ and $~\nu ^{\pm }_{z}$ are adaptive, viz., they are re-evaluated at each time step to ensure that dissipation occurs near the smallest scales of the grid in order to maximize the inertial range. This is necessary because the turbulent amplitudes change by orders of magnitude over the course of the simulations, thus changing the dissipation scale for a given (fixed) hyperviscosity significantly. The method is explained in more detail in Appendix B.

2.2.1 Simulation parameters

In the expanding RMHD equations, there are three ratios of time scales that will prove important for the dynamics. We will define these in more detail below, but feel it useful to introduce the notation here: $\chi _{A}$ will denote the usual ratio of Alfvénic to nonlinear time scales (Goldreich & Sridhar 1995; Mallet, Schekochihin & Chandran 2015); $\chi _{\rm exp}$, the ratio of the expansion to nonlinear time scales; and $\varDelta = \chi _{\rm exp}/\chi _{A}$, the ratio of the expansion to Alfvénic time scales:

(2.11a–c)\begin{equation} \chi_{A} \doteq \dfrac{\tau_{A}}{\tau_{\rm nl}}, \quad \chi_{\rm exp} \doteq \dfrac{\tau_{\rm exp}}{\tau_{\rm nl}}, \quad \varDelta \doteq \dfrac{\tau_{\rm exp}} {\tau_{A}}. \end{equation}

Note that, because $\tau _{A}$ and $\tau _{\rm nl}$ are both proportional to the size of a given structure, they decrease towards smaller scales, while $\tau _{\rm exp}$ is independent of scale. Thus, the effect of expansion decreases in importance towards smaller and smaller scales in the turbulence, eventually reverting to standard RMHD ($\chi _{\rm exp}\rightarrow \infty$, $\varDelta \rightarrow \infty$).

Because of the rescaling symmetry of the RMHD equations, aside from resolution and dissipation properties, two of these three parameters set the important parameters of a given simulation. It is most natural to set $\chi _{\rm exp }$ and $\chi _{A}$ via the initial conditions (discussed below) and fix the ratio of the box-scale Alfvén frequency ($\omega _{A,{\rm box}}= 2{\rm \pi} v_{A}/L_{z}$) to the expansion rate,

(2.12)\begin{equation} \varDelta_{\rm box} \doteq \frac{\omega_{A,{\rm box}}}{\dot{a}/a}= \dfrac{2{\rm \pi}}{L_z}\dfrac{v_{A0}}{\dot{a}}, \end{equation}

where the second expression accounts for the fact that $\varDelta _{\rm box}$ remains constant throughout the evolution because $v_{A}\propto 1/a$.

We choose $\varDelta _{\rm box}$ in the simulation by reference to the conditions observed around the Alfvén radius $R_{A}$ (Kasper et al. 2021). Assuming purely parallel fluctuations on a radial field, taking $\dot {a}/a = U/R$ (see (2.1)),³ and ignoring violations to the Taylor hypothesis so that fluctuations of parallel scale $\ell _z$ yield a spacecraft frequency $f\approx U/\ell _z$, we find that

(2.13)\begin{equation} \varDelta \approx 19.6 \,\mathcal{M}_{A}^{{-}1}\frac{f}{10^{{-}4}\ {\rm Hz}}\frac{R}{18 R_\odot} \left( \frac{U}{400\ {\rm kms}^{{-}1}}\right)^{{-}1}.\end{equation}

Here $\mathcal {M}_{A}=U/v_{A}$ and the normalization to $f=10^{-4}\ {\rm Hz}$ is chosen as the minimum value measured around these radii, which sits well below the measured correlation scales (e.g. Kasper et al. 2021; Zank et al. 2022). Based on this estimate, we set

(2.14)\begin{equation} \varDelta_{\rm box}=10 \end{equation}

for all simulations. Because $\mathcal {M}_{A}\propto R$ and $U\sim {\rm const.}$ in the super-Alfvénic wind, the minimum resolved parallel frequency of ${\simeq }5\times 10^{-5}\ {\rm Hz}$ (i.e. that corresponding to $\varDelta _{\rm box}=10$ from (2.13)) remains constant as the simulation evolves. Given that the correlation scale is observed and expected to increase with $R$, it may be important to understand the effect of lowering $\varDelta _{\rm box}$ in future work.

Note that the choice of $\varDelta _{\rm box}$ can be equivalently understood as setting the resolution in $k_z$ of the simulation: with infinite spatial resolution, a longer box, which contains lower $k_{z}$ modes, is identical to a shorter box with smaller $\varDelta _{\rm box}$. It is also of note that there exist $k_{z}=0$ two-dimensional (2-D) modes, which do not propagate, unlike the other modes in the box. While these are, in some respects, an artefact of the expanding box's periodic boundary conditions, we argue below that they are capturing important physical effects and should not be artificially excluded. More discussion of the choice of $\varDelta _{\rm box}$ and its possible effect on the turbulence is given in § 5.

2.2.2 Initial conditions

Rather than realistically simulate a patch of solar wind as it propagates outwards, the goal of this work is to distill and understand theoretically the key physical features of reflection-driven turbulence. Therefore, our initial conditions are idealized and designed to understand the model itself, on the assumption that this is a prerequisite for understanding the physical processes it attempts to represent. Anticipating the result that the correlation scales of the turbulence will increase significantly as it evolves, it is thus important to start with fluctuations on scales well below the box scale in order to avoid artificially constraining the system's evolution. We choose to obtain the initial $\tilde {\boldsymbol {z}}^+$ fields from a balanced RMHD simulation evolved into its statistically stationary turbulent state, loosely motivated by the idea that outwards Alfvénic fluctuations could ‘escape’ into the corona through an effective high-pass filter from a region of nearly balanced stronger turbulence (van Ballegooijen et al. (2011); although, of course, the EBM is formally valid only outside the Alfvén point by which point the turbulence will have evolved).⁴ The forcing of this balanced simulation is local in Fourier space, acting on all the modes within the ring $k_{\perp } \in 2{\rm \pi} /L_{\perp } [ 99.5, 100.5 ]$ and $\vert k_z \vert = 2{\rm \pi} /L_{z}$, and is designed so as to keep the rate of injection of energy constant with the amplitude needed for critical balance (Goldreich & Sridhar 1995). This creates initial fluctuations with a correlation scale modestly above the forcing scale, with a perpendicular correlation length $L_+\approx L_{\perp }/75$ and parallel correlation length ${\approx }L_{z}$. In the infra-red range (scales larger than the perpendicular correlation length), the initial energy spectrum scales approximately as ${\propto } k_\perp$ in accordance with theoretical expectations (Schekochihin 2022). We use the $\tilde {\boldsymbol {z}}^+$ field thus obtained to initialize the $\tilde {\boldsymbol {z}}^{-}$ one by setting $\tilde {\boldsymbol {z}}^{-}=-\kappa \tilde {\boldsymbol {z}}^+$ with $\kappa$ such that $1-\sigma _c=1\times 10^{-4}$ (this choice is not of great importance because the system rapidly self-adjusts).

Given this choice of a spectrum of fluctuations, the only remaining parameter of interest is the RMHD fluctuation amplitude, which, as shown below, has a strong impact on the turbulence evolution. Based on rescaling symmetry discussed above, this amplitude should be thought of as controlling the ratio of the nonlinear time scale $\tau _{\rm nl}^{\mp }\sim (k_{\perp } z^{\pm })^{-1}=a^{-3/2}(\tilde {k}_{\perp }\tilde {z}^+)^{-1}$ to the linear time scales $(k_{\|}v_{A})^{-1}$ and $a/\dot {a}$, as opposed to directly setting the physical turbulent amplitude $z^+/v_{A}$ (or $|\boldsymbol {B}_{\perp }|/\bar {B}$ or $|\boldsymbol {u}_{\perp }|/v_{A}$). Accordingly, we set

(2.15a,b)\begin{equation} \chi_{\rm exp0} \doteq \frac{k_{\perp0} z^+_{rms0} }{\dot{a}/a} \quad\text{and} \quad\chi_{A0} \doteq \frac{k_{\perp0} z^+_{rms0} }{k_{z0} v_{A}} \end{equation}

as simulation parameters by rescaling $\tilde {\boldsymbol {z}}^+$ by the required amount. Here $k_{\perp 0}$ and $k_{z0}$ are the initial inverse correlation lengths, $z^+_{rms0}$ is the initial root-mean-square (r.m.s.) fluctuation amplitude and the ratio $\chi _{\rm exp0}/\chi _{A0}=\varDelta _{\rm box}$ is fixed to be 10 for all simulations as described above (i.e. rescaling $\tilde {\boldsymbol {z}}^+$ sets both $\chi _{A0}$ and $\chi _{\rm exp0}$ because we have already fixed $\varDelta _{\rm box}$). We shall see that because they have stronger nonlinearity, simulations with larger $\chi _{\rm exp 0}$ remain in the strongly nonlinear regime for longer, thus displaying more clearly the relevant power-law behaviour and clarifying the analysis. Most figures and discussion will thus focus on the highest-$\chi$ case run, which has $\chi _{\rm exp0}=960$ ($\chi _{A0}=96$) and a resolution $n_\perp ^2\times n_z=1536^2\times 256$. This value of $\chi _{\rm exp0}$ is rather large compared with the solar wind around $R_{A}$ at the correlation scale of the turbulence (in § 7 we estimate that $\chi _{\rm exp0}\approx 60$ near the Alfvén point in the conditions observed by Kasper et al. 2021); however, $\chi _{\rm exp0}$ likely varies significantly between streams. We have run a series of simulations down to $\chi _{\rm exp0}=0.75$ and will show some of these for comparison. Note that various previous works (e.g. Dong et al. 2014; Montagud-Camps, Grappin & Verdini 2020; Grappin et al. 2022) have denoted $\chi _{\rm exp}$ by $\epsilon ^{-1}$, with $k_{\perp 0}$ taken to be $2{\rm \pi} /L_\perp$, exploring expanding turbulence for $\epsilon \gtrsim 0.2$ ($\chi _{\rm exp}\lesssim 5$).

This method of constructing initial conditions, while straightforward and well controlled, has the downside of placing the plasma into an artificial ‘super-critically balanced’ state ($\chi _{A}\sim k_{\perp }z^+/k_{\|}v_{A}>1$). The consequence is that, over a relatively short transient initial phase as the fields start evolving nonlinearly, neighbouring planes along the $\hat {\boldsymbol {z}}$ direction decorrelate and develop small parallel scale fluctuations until $\chi _{A}\sim 1$, establishing critical balance. This transient process generates a flat $k_\|$ spectrum (white noise) up to the parallel scale at which $k_\| v_{A}$ balances the nonlinear mixing, which is the Fourier-space hallmark of critical balance (Schekochihin 2022). This process occurs over a time scale comparable to the nonlinear time at each scale, which is rapid compared with the time it takes the system to decay and change regimes, so we believe this choice does not strongly impact our results. However, future work should explore the effect of this choice, other initial conditions and $\varDelta _{\rm box}$ in more detail in order better understand the impact of our choices.

3 Basic evolution

Starting from the initial conditions described above, we evolve the system with $a$ (equivalently, with time) up to $a=1000$. While this would correspond, in principle, to an extremely large physical radius [$R\approx 1000R_{A}\approx 84\ {\rm au} (R_{A}/18 R_{\odot })$], we reiterate that we are deliberately exploring more extreme parameters in order to better characterize the physics of reflection-driven turbulence. For more realistic initial conditions with lower $\chi _{\rm exp0}$, the behaviour and transitions we describe below will occur at smaller $a$.

As illustrated in figure 1, which shows important aspects of how simulations with different initial amplitudes ($\chi _{A0}$) evolve with $a$, the system's evolution is naturally divided into two distinct phases, discussed separately in §§ 4 and 5 below. Following a short initial transient, when $\boldsymbol {z}^{-}$ and the parallel scales rapidly adjust (see above), the first ‘imbalanced’ phase involves turbulence where the normalized cross-helicity, or imbalance,

(3.1)\begin{equation} \sigma_c \doteq \dfrac{\tilde{E}^+-\tilde{E}^{-}}{\tilde{E}^++\tilde{E}^{-}} = \frac{\tilde{E}^{c}}{\tilde{E}} \end{equation}

is almost maximal (unity), as in the initial conditions. In the strong nonlinear regime ($\chi _{A0}=96$; solid lines in figure 1a), the turbulent energy decays as $\tilde {E}\approx \tilde {E}^+\propto a^{-1}$, signalling turbulent heating of the plasma. In contrast, in the second ‘magnetically dominated’ or balanced phase, which starts at around $a\approx 80$ for the $\chi _{A0}=96$ simulation, $\sigma _{c}$ approaches $0$ with $\tilde {E}^+\approx \tilde {E}^{-}$, and surprisingly, $\tilde {E}$ starts growing in time. This is a consequence of the system developing a large negative normalized residual energy,

(3.2)\begin{equation} \sigma_r\doteq\dfrac{\tilde{E}^{u}-\tilde{E}^{b}}{\tilde{E}^{u}+\tilde{E}^{b}} = \frac{\tilde{E}^{r}}{\tilde{E}}, \end{equation}

which, as seen from (2.8), can cause $\tilde {E}$ to grow as $\tilde {E}\propto a$ (as observed) in the absence of dissipation. We show this evolution graphically with the ‘circle plot’ in figure 1(b). This illustrates the evolution of $\sigma _c$ and $\sigma _r$ during the radial transport (Bruno et al. 2007), which are constrained by the relationship between $\tilde {E}^\pm$, $\tilde {E}^u$ and $\tilde {E}^b$ to lie within the circle $\sigma _c^2+\sigma _r^2=1$. The fact that the evolution remains near the edge of the circle indicates that the fields maintain a high level of ‘Elsässer alignment’ between $\tilde {\boldsymbol {z}}^+$ and $\tilde {\boldsymbol {z}}^{-}$, with

(3.3)\begin{equation} \sigma_{\theta}\doteq \dfrac{\langle \tilde{\boldsymbol{z}}^+{\cdot}\tilde{\boldsymbol{z}}^{-} \rangle}{\langle |\tilde{\boldsymbol{z}}^+|^{2}\rangle^{1/2}\langle|\tilde{\boldsymbol{z}}^{-}|^{2}\rangle^{1/2}}=\dfrac{\sigma_{r}} {\sqrt{1-\sigma_c^{2}}} \end{equation}

close to $-1$ in the later stages of the simulation (the isocontours of $\sigma _{\theta }$ are shown by solid lines in figure 1). This strong alignment is likely primarily a consequence of the reflections, which generate $\tilde {\boldsymbol {z}}^{-}$ fluctuations that are perfectly aligned with $-\tilde {\boldsymbol {z}}^+$, although the mutual shearing of the Elsässer fields is also known to generate aligned fluctuations even in homogeneous Alfvénic turbulence (Chandran et al. 2015). The simulation's evolution bears a striking resemblance to the joint distribution of normalized cross-helicity and residual energy observed in highly Alfvénic fast-solar-wind streams (Bruno et al. 2007; Wicks et al. 2013; D'Amicis et al. 2021), providing good evidence that, despite the drastic approximations involved with our model, it captures some of the key physics of solar-wind turbulence.

Figure 1. (a) Radial evolution of wave-action energies $\tilde {E}^+$ (red lines) and $\tilde {E}^-$ (blue lines) for three simulations with different amplitude initial conditions. Solid lines show our highest-amplitude $\chi _{\rm exp0}=960$ ($\chi _{A0}=96$) simulation, dash-dotted lines show the $\chi _{\rm exp0}=7.5$ ($\chi _{A0}=0.75$) simulation and dotted lines show the $\chi _{\rm exp0}=0.75$ ($\chi _{A0}=0.075$) simulation in the weak regime. We normalize each curve to its initial $\tilde {E}^+$ to facilitate comparison and the dotted-grey lines indicate various power laws for reference (see text). (b) Parametric representation of $\sigma _r$ and $\sigma _c$ during the evolution of the $\chi _{A0}=96$ simulation. The colours (on a logarithmic scale) indicate the normalized radial distance $a$. Solid lines represent contours of constant $\sigma _\theta$ as labelled (see text).

The properties of the turbulence change dramatically between the two phases, as illustrated by the perpendicular snapshot of $\tilde {\boldsymbol {z}}^{\pm }$ shown in figure 2. Most obviously, the turbulence dramatically increases in scale with time, starting from the very small scales of the initial conditions (a,b) to reach nearly the box scales by the latest times (e,f). We will argue below that this is a consequence of the anomalous turbulent growth of ‘wave-action anastrophy’ during the imbalanced phase, which significantly constrains the turbulence as it decays, forcing it to rush to larger scales and form a split cascade. At early times, the structures in $\tilde {\boldsymbol {z}}^+$ and $\tilde {\boldsymbol {z}}^{-}$ are rather different, with different dominant scales, but as the turbulence enters the magnetically dominated phase (c,d) the two become more similar as it becomes balanced. A key change (not shown in figure 2) is that the turbulence becomes more two dimensional at larger $a$, with structures across a wide range of $k_{z}$ at earlier times (a,b) giving way to predominantly $k_{z}=0$ modes by the $a=250$ snapshot shown in (e,f). While true $k_{z}=0$ modes are of course an artefact of the periodicity of the EBM, their key feature as pertains to reflection turbulence is that they are expansion dominated and do not propagate, unlike Alfvén waves. Since this is the case for any sufficiently long-wavelength mode, even in non-periodic settings or the real solar wind (specifically, those with $\varDelta = k_{z}v_{A}/(\dot {a}/a)<1/2$; see § 5), we argue that these dynamics are physical and likely have already been observed in the solar wind (Tu & Marsch 1991). As seen also in figure 1(a), there is little turbulent heating in this phase, which (we will show) occurs because the circular structures approach local nonlinear ‘Alfvén vortex’ solutions (Petviashvili & Pokhotelov 1992; Alexandrova 2008), which slows down their evolution significantly, impeding their dissipation.

Figure 2. Snapshots of the Elsässer fields $|\tilde {\boldsymbol {z}}^+|$ (a,c,e) and $|\tilde {\boldsymbol {z}}^{-}|$ (b,d,f) across the full box in the plane perpendicular to the mean magnetic field for three different radial distances. Panels (a,b) illustrate $a=5$ during the imbalanced decay phase; (c,d) show $a=50$, which is shortly before the transition to the balanced phase; (e,f) show $a=250$ in the balanced, magnetically dominated regime. This simulation has a resolution of $n_\perp ^2\times n_z=8192^2\times 256$ and is initialized by progressively refining the $n_\perp ^2 \times n_z= 1536^2 \times 256$ simulation that was run from $a=1$ to $a=1000$.

We now explore the two phases in more detail, attempting to diagnose and understand key features of their turbulence to make detailed predictions for solar-wind observations.

4 Imbalanced phase

In this section we explore the turbulence in the imbalanced phase of the simulations, which applies when $\tilde {z}^+\gg \tilde {z}^{-}$, for $a\lesssim 50$ in the $\chi _{A0}=96$ simulation (see figure 1). Based on figures 1 and 2, the key features of this phase that we wish to understand are (i) the power-law evolution of $\tilde {E}^{\pm }$, which sets the overall heating (turbulent-decay) rate as a function of radial distance; and (ii) the cause of the significant increase in the fluctuations’ scale during their evolution. To interpret the basic time evolution of $\tilde {E}^+$ and $\tilde {E}^{-}$, we first (§ 4.1) review and assess phenomenological ideas based on Dmitruk et al. (2002), which have been used in a number of previous works to predict and understand reflection-driven turbulence both inside and outside the Alfvén point (Verdini & Velli 2007; Chandran & Hollweg 2009; Verdini et al. 2010; Chandran & Perez 2019). While the phenomenology is consistent with some general features of the observed time evolution (§ 4.1) and spectra (§ 4.2), we will find some important differences that we cannot, at this point, satisfactorily explain. Whether these signal fundamental issues with the theoretical basis of the model, or just more minor discrepancies, remains unclear. In this discussion we will see that feature (ii) (the rapid increase in the scale of the fluctuations) happens to not influence the decay, so it can be discussed separately. We argue in § 4.4 that this feature arises from the surprising property of anomalous turbulent wave-action anastrophy growth, which constrains the dynamics and forces $\tilde {z}^+$ to rush to large scales as it decays via a split cascade.

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. 1989; Lithwick, Goldreich & Sridhar 2007; Chandran & Perez 2019; Schekochihin 2022), 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. 2002):

(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 2007; Chandran & Hollweg 2009; Verdini et al. 2010) 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 1941), 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.1.1 Numerical results

Consider first the figure 1(a), focusing on the decay (growth) of $\tilde {E}^+$ ($\tilde {E}^{-}$) for the $\chi _{A0}=96$ ($\chi _{\rm exp0}=960$) simulation (solid lines), which undergoes a long period of power-law evolution before reaching the balanced regime. We see that $\tilde {E}^{-}\propto a^{2}$, which is significantly faster than the simplest prediction from (4.3) with $\tilde {\lambda }_+\sigma _{\theta }\sim {\rm const.}$ (yielding $\tilde {E}^{-}\propto a^{1}$). While this is perhaps not surprising, since, as seen in figure 2, the fluctuations’ scales are increasing rapidly with time (thus, presumably increasing $\tilde {\lambda }_+$), we have not identified a clear candidate for providing the additional factor of $a^{1/2}$ in (4.3).⁵ The $\tilde {E}^+$ decay, in contrast, matches the simplest prediction of (4.4), with $\sigma _{\theta }\tilde {\lambda }_+/\tilde {\lambda }_{-}\approx 1$ and $\tilde {E}^+\propto a^{-1}$. This feature seems robust across different initial conditions with sufficiently high $\chi _{\rm exp0}$ and suggests physically that the dominant scale of $\tilde {z}^{-}$ that advects $\tilde {z}^+$ to cause dissipation ($\tilde {\lambda }_{-}$) is the same as that which governs the evolution of $\tilde {z}^{-}$ ($\tilde {\lambda }_+$). The reason for such a correspondence is not immediately obvious but may be related to the fact that the scales that control the growth of $\tilde {z}^{-}$ are also coherent with $\tilde {z}^+$ (being driven by reflection), and thus, most effective at advecting and dissipating $\tilde {z}^+$. Another, non-exclusive possibility is that the $\tilde {z}^{-}$ that is most effective at advecting $\tilde {z}^+$ has a ${\sim }k_{\perp }^{-3}$ spectrum (as motivated above), which would yield a nonlinear turnover time ($\tau _+\sim a^{3/2}\tilde {\lambda }_{-}/\tilde {z}^{-}$) that is independent of scale. Perhaps also of note is that a self-similar power-law decay is possible in (4.4) only if $\sigma _{\theta }\tilde {\lambda }_+/\tilde {\lambda }_{-}$ is constant.

The lower-amplitude simulations, with $\chi _{A0}=0.75$ ($\chi _{\rm exp0}=7.5$; dash-dotted lines) and $\chi _{A0}=0.075$ ($\chi _{\rm exp0}=0.75$; coloured dotted lines) behave rather differently. The $\chi _{A0}=0.75$ shows a small amount of decay in $\tilde {E}^+$, while the $\chi _{A0}=0.075$ case shows almost none, and $\tilde {E}^{-}$ grows much more rapidly and is not a power law in either simulation. We will show in § 5 that this behaviour is effectively just the linear growth of long-wavelength $\tilde {z}^{-}$ modes, which are $k_{z}=0$ modes seeded from the initial conditions in the simulation. The linear growth of such modes is significantly faster than the nonlinear prediction (4.3), so the system reaches the balanced regime at smaller $a$ (equivalently, the nonlinear prediction is $\tilde {z}^{-}\sim \tilde {z}^+/\chi _{\rm exp}$ and $\chi _{\rm exp}$ is not large initially). The lack of $\tilde {E}^+$ decay is a consequence of the turbulence being weak, or, in the case of the $\chi _{A0}=0.75$ simulation, rapidly becoming so, because $\chi _{A}\sim (\tau _{-})^{-1}/(k_{\|}v_{A})\propto a^{-1/2}$ for fixed $\tilde {z}^+$ and $k_{\|}$. We have observed generically that weak turbulence in the EBM exhibits almost no nonlinear decay, behaving effectively as a collection of linear modes. However, we caution that key aspects of the expanding-box approximation are not valid for modes in the weak regime, and its predictions for how $\tilde {\boldsymbol {z}}^{-}$ is forced via randomly phased $\tilde {\boldsymbol {z}}^+$ are likely incorrect.⁶ Further work is needed to understand these issues, but weak-turbulence EBM results should be treated with caution.

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. ⁷ 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. (1989), Verdini, Velli & Buchlin (2009), Perez & Chandran (2013) and Chandran & Perez (2019). 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. 2009). 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. (1989) and extended to more realistic anisotropic turbulence by Perez & Chandran (2013), 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. 2009; Perez & Chandran 2013; Chandran & Perez 2019; Squire, Chandran & Meyrand 2020) 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 2016; Lugones et al. 2019), 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 prediction⁸ $\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 2019). 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. 1989; Lithwick et al. 2007; Perez & Chandran 2013; Chandran & Perez 2019; Schekochihin 2022), 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. (2019) and Yang et al. (2023) 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.⁹ 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. 2009). 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.¹⁰ 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 }^{-}$.¹¹ 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 1941). 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 1978; Alexakis & Biferale 2018). 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 1984; Biskamp & Welter 1989). 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 1995):

(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 2017) and observed in high cross-helicity solar-wind streams (Smith et al. 2009; Coburn et al. 2015).

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 Balanced, magnetically dominated phase

With $\tilde {z}^+$ decaying while $\tilde {z}^{-}$ grows, it is clear that the imbalanced phase must inevitably end as the fluctuations approach the balanced regime with $\tilde {z}^{-}\sim \tilde {z}^+$. Indeed, recall that the phenomenology of § 4.1 predicted $\tilde {z}^{-}\sim \tilde {z}^+/\chi _{\rm exp}$, where $\chi _{\rm exp}$ is the ratio of expansion to nonlinear times (4.7), which is necessarily a decreasing function of time. Thus, $\chi _{\rm exp}\sim 1$ marks the end of the imbalanced phase. In figure 1 we saw that the wave-action energy starts growing in time, with $\tilde {E}\propto a$, characterized by magnetically dominated fluctuations ($\sigma _{r}<0$), and with minimal turbulent dissipation into heat. It is the purpose of this section to understand the important properties of this balanced, magnetically dominated phase, making predictions for in-situ observations at large distances from the sun.

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 1980; Hollweg 1990; Hollweg & Isenberg 2007); 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 1980). 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 1991; Grappin et al. 1993). 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 1980; Hollweg & Isenberg 2007), 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. 2015) – a vivid example of self-organization (Horiuchi & Sato 1985). They can be understood using a classical variational argument (Matthaeus & Montgomery 1980). 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 1992; Alexandrova 2008; Lion, Alexandrova & Zaslavsky 2016). 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. 2019). 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 1991). 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. 2019), 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 1984; Zhou et al. 2019). 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 1954; Schuurman, Bobeldijk & de Vries 1969). 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 Relationship to previous literature

Naturally, the subject of solar-wind turbulence and transport has a history as long as the space-age itself, and there exists a wide literature attempting to model, and ultimately predict, observations of turbulence, heating, wind speed and magnetic fields. Given this history, in this section we feel it is useful to highlight the similarities and differences between our work and previous literature. We mostly restrict the discussion to super-Alfvénic solar wind, focusing only on how our model and results relate to previous works, rather than providing a comprehensive literature review. The literature can be divided into two broad groups: first, discussed in § 6.1, those models where turbulent nonlinearities are approximated using a one-point closure; and second (§ 6.2), where the turbulence is captured directly via fluid equations. The majority of our results fall into the second category, although our analysis and discussion is necessarily informed by the first.

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 (1989), Zhou & Matthaeus (1989), Velli et al. (1989) followed by Matthaeus et al. (1994) and Matthaeus, Zank & Oughton (1996). 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. 2014; Usmanov et al. 2018).

6.1.1 Shear-driven models

The idea that large-scale inhomogeneities in the solar wind would provide free energy from unstable velocity shear to drive turbulence continuously as the wind expands was first proposed by Coleman (1968). This idea later gained considerable momentum when two-scale transport theories for solar-wind fluctuations were developed. Today it is a leading paradigm for describing super-Alfénic solar-wind evolution. With the addition of interstellar pickup ions, which provide another source of energy injection into the distant solar wind, these models found good agreement with spacecraft data for reasonable choices of parameters (Zank et al. 1996; Matthaeus et al. 1999), explaining their subsequent popularity (Breech et al. 2008). All of these models include so-called ‘mixing’, which corresponds to the exchange of energy between Elsässer fields through expansion and compression, as well as large-scale sheared flows. These models have been therefore referred to as MECS. Note that most incarnations of such models ignore Alfvénic dynamics, a distinguishing feature from the reflection-driven models discussed below. In recent years, more complete models based on similar ideas include solutions in full 3-D geometry (e.g. Usmanov et al. 2018; Chhiber et al. 2021) and/or more detailed turbulence closures to model different types of fluctuations (Zank et al. 2011; Adhikari et al. 2015; Zank et al. 2017).

A main conclusion of this body of literature is that spherical expansion alone, and therefore reflections, cannot explain turbulent transport – shear is needed. It is interesting to revisit this conclusion in view of the results discussed in § 5. Indeed, in our reflection-driven simulations, as the system transitions into the balanced phase and becomes dominated by slowly merging Alfvén vortices, the magnetic energy follows a radial evolution close to the one predicted by shear-driven models (Zank et al. (1996, 2017); see § 7.2.6 below for further details). How is it that two different models, containing different physical ingredients, give similar results? The answer lies in the assumptions of particular turbulence symmetries that underpin the MECS-based turbulent transport models. Chiefly, the assumption that the normalized residual energy $\sigma _r$ is a constant and that magnetic tension, which couples velocity and magnetic field, can be ignored. These two assumptions are inconsistent with each other and may lead to non-physical conclusions. Indeed, because of the density stratification, in the absence of Alfvénic coupling, magnetic and kinetic energy cannot decay at the same rate in an expanding medium. This explains the development of highly magnetically dominated states that impede the nonlinear decay through the formation of Alfvén vortices, as we observe in our simulation (see § 5.3). Because this impeding effect is absent by construction in the MECS transport models, energy replenishment is needed to compensate for the assumed strong turbulent damping, leading to the conclusion that shear and pickup ions are required. Our results challenge this interpretation. We argue, supported by direct numerical simulations, that although large-scale shear and pickup ions are certainly present in the solar wind, they are not necessary to explain the evolution of magnetic fluctuation amplitudes at large radial distances. It should be noted that more recent versions of similar two-scale closures do include both Alfvénic physics and a transport equation for the residual energy (Zank et al. 2011). In this model, however, the injection of reduced energy is linked to shear and is modelled by a free parameter that is adjusted to fit observations (Adhikari et al. 2015; Zank et al. 2017, 2021). These differences in the treatments of residual energy account for the different results of our models.

6.1.2 Reflection-driven models

The key ideas required for a reflection-driven turbulence phenomenology were first discussed in Velli et al. (1989). Reflection-driven two-scale transport models (Matthaeus et al. 1999; Cranmer 2005; Suzuki & Inutsuka 2005) include spherical expansion but ignore compression and shear, but – unlike most MECS-based models – include Alfvénic physics (in this sense, they might be termed AMEN models – ‘Alfvénic, mixing, expansion and nothing else’). The phenomenological model we have used to describe our numerical results falls into this AMEN category. Most previous reflection-driven models have been mainly developed to describe the physics of the lower corona, inside the Alfvén critical point, where sub-Alfvénic physics like super-radial expansion are very important for a reasonable description of the turbulent transport. Our model excludes such physics and also differs from previous ones in that it includes the radial evolution of the turbulence correlation scales (historically ignored in reflection-driven models, but kept in MECS-based models) and the existence of different regimes that shut off the turbulent cascade at larger radii once $\chi _{\rm exp}\sim 1$. It also emphasizes the crucial role played by the residual energy in the approach to the nonlinear Alfvén vortex solutions.

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. 2002; Matsumoto & Suzuki 2012; Perez & Chandran 2013; Shoda et al. 2019). 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. 1993; Dong et al. 2014). 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 (2013) and van Ballegooijen & Asgari-Targhi (2016). 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. (2009, 2012) 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. 2020; Johnston et al. 2022). Another important difference is that most previous works (e.g. Dong et al. 2014; Verdini & Grappin 2015, 2016; Montagud-Camps et al. 2018) 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 2013; Wicks et al. 2013). 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. (2020) and Grappin et al. (2022) 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. 2020). 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. (2020), Montagud-Camps et al. (2020) and Grappin et al. (2022), 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 Solar-wind observations

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. 2021). 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. 2021), 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. 1987; Tu & Marsch 1995). 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 1980; Grappin et al. 1982; Chen et al. 2011; Schekochihin 2022). Some suggestions invoke interactions between different streams as the dominant influence (Roberts et al. 1992; Matthaeus et al. 2004; Adhikari et al. 2015), or parametric decay of Alfvén waves (Goldstein 1978), 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. 1989; Zhou & Matthaeus 1990; Oughton & Matthaeus 1995; Verdini et al. 2009; Grappin et al. 2022), 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. 2002; Verdini & Velli 2007) 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. 2007; Chandran 2021), in which slower streams arise because more of the outward-fluctuation energy is dissipated at low altitudes (Halekas et al. 2023), 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. 2021), 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 1998; Bruno et al. 2007; Bruno & Carbone 2013; Wicks et al. 2013). 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. 1998; Bruno et al. 2007), 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 1990; Oughton & Matthaeus 1995; Hollweg & Isenberg 2007; Grappin et al. 2022), 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.3 Fluctuation spectra and the $1/f$ range

Many years of observations have shown that magnetic fluctuations in the solar wind display a $1/f$ slope at low frequencies, differing from the steeper $f^{-3/2}$ or $f^{-5/3}$ scalings observed at higher frequencies in the inertial range (Goldstein, Burlaga & Matthaeus 1984; Matthaeus & Goldstein 1986; Chen et al. 2020). There is currently no consensus on the origin of this $1/f$ range. Suggestions range from its origin in the low corona (Matthaeus & Goldstein 1986), implying that it is the energy reservoir that feeds the solar-wind turbulent cascade (Matthaeus et al. 1994), to it being the result of spherically polarized fluctuations growing to amplitudes larger than one (Matteini et al. 2018) or parametric decay of compressive fluctuations (Chandran 2018; Davis et al. 2023; Huang et al. 2023). Numerous studies starting from Velli et al. (1989) have also shown that reflection-driven turbulence can naturally create $1/f$ spectra in both the parallel (Verdini et al. 2012) and perpendicular (Perez & Chandran 2013; Chandran & Perez 2019) directions. Our results agree with the latter¹² through the mechanism described in Velli et al. (1989), Perez & Chandran (2013) (see § 4.2). In line with previous works (Perez & Chandran 2013; Chandran & Perez 2019), we find that the spectral scaling of $\tilde {z}^{-}$ is steeper than that of $\tilde {z}^+$ through this range, scaling as $\tilde {E}^{-}\propto k_{\perp }^{-3/2}$ in our simulations; this is similar (though not identical) to that observed in situ (Tu, Marsch & Rosenbauer 1990; Wicks et al. 2013), although this general signature is not unique to the reflection-turbulence model (Chandran 2018; Matteini et al. 2018). In addition, since the $1/f$ range in the model is generated by the turbulence during the imbalanced phase, at similar radii, we would expect $\chi _{\rm exp}\gtrsim 1$ regions to exhibit a wider $1/f$ range than $\chi _{\rm exp}\lesssim 1$ regions. If we further apply the hypothesis discussed above, that fast-wind streams have higher $\chi _{\rm exp0}$ than slow-wind streams, this would naturally explain the well-known observation that the size of the $1/f$ range correlates with wind speed (Tu, Marsch & Thieme 1989; Tu et al. 1990). In this context, it is also worth clarifying that the extremely wide $k_{\perp }^{-1}$ range seen in figure 3 is again a consequence of the extreme parameters of the simulation (its small initial scale and long imbalanced phase). Finally, the general ideas naturally explain the results of Wicks et al. (2013) that in those (rare) regions with $\sigma _r\simeq \sigma _c\sim 0$, there is no significant $1/f$ range, since (given figure 1) such regions are presumably strongly influenced by physics that is unrelated to reflection-driven turbulence.

7.2.4 Inverse energy transfer and the split cascade

The most significant qualitative difference between our energy spectra and previous results is the inverse energy transfer of $\tilde {E}^+$ caused by anomalous growth of wave-action anastrophy. This forces the decay to proceed via a split cascade, shifting the $\tilde {E}^+\propto k_{\perp }^{-1}$ range to larger scales with time in the co-moving frame as it grows out of a positive-slope infrared spectrum at yet larger scales. The feature is interesting in light of recent observations showing that the $1/f$ spectrum does not extend to the largest available scales, especially near the sun (Kasper et al. 2021; Huang et al. 2023), instead developing as the wind propagates outwards (Davis et al. 2023). The surprising, non-trivial prediction of our model is that the correlation scales of the fluctuations, which lie towards the large-scale side of the $1/f$ range, should increase with $R$ faster than expansion (i.e. increase in the co-moving frame). In addition, the split cascade itself may be directly observable, with the prediction that the cascade of $z^+$ should switch from forward to inverse at the largest scales in imbalanced turbulence (see figure 6). Interestingly, back transfer of energy from small to large scales has been observed in $z^+$ in highly imbalanced streams (Smith et al. 2009; Coburn et al. 2015), although since these observations seem to pertain to smaller scales (where we observe a forward cascade of both $z^+$ and $z^{-}$; figure 6), they may be unrelated.

This inverse energy transfer could have broader implications for solar-wind turbulence and acceleration, particularly if similar physics also applies in sub-Alfvénic regions. Close to the Sun, the large gradient of the Alfvén speed around the transition region should prevent low-frequency Alfvén waves launched from the chromosphere from propagating outwards to large radial distances (Hollweg 1978; Leroy 1981; Velli et al. 1991; Réville et al. 2018). If the chromospheric fluctuations are turbulent and critically balanced, with little power in modes with $v_{A} k_\|>z^{\pm } k_\perp$ (Schekochihin 2022), this high-pass filter would also have the effect of filtering out large scales in the perpendicular direction, leading to small correlation scales at the coronal base. The fact that low-frequency waves end up dominating the solar-wind spectrum is therefore highly non-trivial and naturally suggests that some form of inverse energy transfer is needed to explain the existence of large-scale fluctuations at all. The anomalous growth of wave-action anastrophy could provide one such mechanism.

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. 1982; Freeman 1988; Tu 1988), presumably by the dissipation of fluctuations (Verma, Roberts & Goldstein 1995; Vasquez et al. 2007; Chen et al. 2020; Halekas et al. 2023). Such turbulent heating appears to be less important in slower-wind streams, although results remain controversial (Freeman 1988; Totten, Freeman & Arya 1995; Hellinger et al. 2011). 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. 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 1988; Totten et al. 1995; Hellinger et al. 2011) should signal the importance of other physics, or not. Indeed, more complex models based on a similar phenomenology (Cranmer 2009; Chandran et al. 2011) reproduce observed temperature profiles reasonably well out to ${\simeq }1\ {\rm AU}$, and Montagud-Camps et al. (2020) 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. 2023).

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. 1994; Zank et al. 2018), stream interactions across a range of radii (e.g. Roberts et al. 1992; Breech et al. 2008) and parametric decay in highly imbalanced regions (Chandran et al. 2018; Shoda et al. 2019).

7.2.6 Alfvén vortices

The final phases of our simulations are characterized by isolated magnetically dominated nonlinear solutions (Alfvén vortices; Alexandrova 2008; Petviashvili & Pokhotelov 2016), in which magnetic tension balances the total pressure. These structures are dominated by expansion, so not expected to be critically balanced in the usual sense (in our simulations they are truly two dimensional), with sharp boundaries and current rings that separate them from the surrounding more quiescent plasma (see figure 8). The validity of our RMHD model is questionable for this phase, since the vortices would presumably involve significant thermal pressure and density perturbations, which might modify their equilibrium and stability. Furthermore, at large radial distances most streams in the ecliptic have a large Parker spiral angle, invalidating our assumption of a mean radial field, although this may be less of an issue at higher latitudes. Despite these issues, we suggest that our results provide a natural explanation for the so-called ‘magnetic field directional turnings’ (MFDTs) observed in low cross-helicity streams (Tu & Marsch 1991), whose origin has challenged a clear theoretical explanation so far (Bruno et al. 2007; Bruno & Carbone 2013). The observed structures are highly magnetically dominated, with an approximate balance between thermal and magnetic pressure and very sharp boundaries in $\boldsymbol {B}$ (Tu & Marsch 1991), just as observed in figure 8. This explanation suggests that MFDTs and Alfvén vortices are the same physical entity whose origin is reflection-driven turbulence. It could be tested by several means such as: (i) a direct fit to observed structures at large radii to (5.7), perhaps with a focus on high-latitude regions where the Parker spiral is less dominant; (ii) by examining their parallel scales, which should satisfy $\varDelta \lesssim 1/2$, thus indicating they are expansion dominated; and (iii) by examining the radial dependence of $\delta B_{\perp }/\bar {B}$, which should grow ${\propto }R$ until reaching large amplitudes ($\delta B_{\perp }/\bar {B}\sim 1$, where our RMHD approximation is no longer valid).

Another interesting connection with past observations concerns the model's predictions for turbulence fluctuation amplitudes as a function of radius, which have been observed by Voyager and Pioneer out to large radial distances. As discussed in § 6, the slow decay of $\delta B/B_0$ compared with WKB predictions has been argued in a range of literature (e.g. Zank et al. 1996; Matthaeus et al. 1999; Breech et al. 2008) to be evidence for the importance of stream interactions and small-scale pickup-ion driven turbulence. However, as seen in figure 1, as a consequence of reaching the low-frequency magnetically dominated, Alfvén vortex solutions, the fluctuation energy decays significantly less rapidly than the WKB prediction, with $\delta B_\perp /\bar {B}$ growing in time. Indeed, the measured evolution of our simulations, with $\delta B_\perp ^2\propto R^{-2.2}$, provides an excellent match to the observed fluctuation amplitude data out to $40\ {\rm AU}$ shown in Zank et al. (1996, figure 4), and fit therein using a model that invoked turbulent driving by stream interactions and pickup ions, which, at large heliocentric distances, gives $\delta B_{\perp }^2\propto R^{-2.3}$ (see also § 6; Zank et al. 2017). While further work with a Parker spiral and compressible effects is clearly needed, we see promising evidence that most aspects of the observed magnetic and velocity fluctuations can be explained without invoking anything other than expansion and Alfvénic (RMHD) physics.

8 Conclusion

This work has presented a detailed computational and phenomenological study of reflection-driven turbulence, which is thought to play a key role in the heating and acceleration of the solar wind (Cranmer 2009), as well as in other magnetized, highly stratified environments such as accretion disk coronae (Chandran et al. 2018). We have approached the problem from the simplest standpoint possible, using the RMHD EBM, which captures Alfvénic (incompressible and perpendicular) dynamics and assumes a constant wind speed $U$ that is faster than the Alfvén speed $v_{A}$. By enabling very high simulation resolutions and clarifying the analysis, this has helped to reveal a rich and non-trivial dynamics that displays features reminiscent of both forced and decaying turbulence paradigms. In order to explore these features in depth, our study has differed from most previous works by deliberately not attempting to match solar-wind parameters, instead focusing on understanding the physical processes from the simplest standpoint possible. While highly idealized, our results can plausibly explain a range of disparate observations from in-situ spacecraft (see § 7), giving us some confidence in the value of the computational approach and the utility of the theoretical framework.

Our most surprising novel results relate to the existence of strong inverse energy transfer, with the decay of the dominant outward-propagating fluctuations proceeding via a split cascade that transfers energy to small and large scales simultaneously. We argue that this results from an anomalous conservation law of the wave-action anastrophy $\tilde {\mathcal {A}}$ (the box-averaged parallel vector potential squared), which can grow due to the effects of expansion in the strongly turbulent system. We provide a heuristic theoretical argument justifying this (§ 4.4) based on linear-wave dynamics and the observation that the Elsässer fields $\boldsymbol {z}^{\pm }$ remain nearly aligned ($\boldsymbol {z}^{-}\propto -\boldsymbol {z}^+$) and anomalously coherent (effectively propagating in the same direction). This latter property, which results from the suppression of collisions between Alfvénic fluctuations, as diagnosed in the simulation via frequency spectra (§ 4.3), leads to a turbulent decay that remains strong even though the minority fluctuations ($\boldsymbol {z}^{-}$) have very low amplitude (Verdini et al. 2009; Perez & Chandran 2013). Using these core ideas, the radial evolution of the energy, imbalance (normalized cross-helicity) and residual energy are analysed via a heuristic phenomenology based on previous works (Dmitruk et al. 2002; Verdini & Velli 2007; Chandran & Hollweg 2009), extended to account for the radial evolution of the different scales of $\boldsymbol {z}^{\pm }$ (§ 4.1). We argue that a key parameter is the ratio of the expansion/reflection time scale to the nonlinear time scale $\chi _{\rm exp}$, which naturally controls the reflection-driven turbulent decay (Dmitruk et al. 2003). Overall the phenomenology provides a reasonable match to most simulation results, although there remain some unresolved discrepancies.

Another result of our work concerns the long-term evolution of the system at large radii, as relevant to the outer heliosphere and regions of slower wind (see below). Our simulations show that super-Alfvénic reflection-driven turbulence is characterized by two distinct phases, separated by the radius at which $\chi _{\rm exp}\sim 1$, which is also where the system becomes balanced ($z^{-}\sim z^+$) and dominated by long-parallel-wavelength modes for which expansion overwhelms the Alfvénic restoring force. From this radius onwards, nonlinear interactions slow significantly as the system cellularizes into a collection of nonlinear Alfvén vortex solutions separated by sharp current-ring boundaries. The structures, which are strongly magnetically dominated, slowly move and merge while their normalized amplitude $|\boldsymbol {B}_\perp |/|\bar {\boldsymbol {B}}|$ grows significantly faster than the usual WKB prediction due to expansion.

8.1 Observations

Despite the simplicity of our RMHD expanding box and the many important physical effects that are unjustifiably neglected (see § 8.2 below), its predictions seem to explain a range of different well-known solar-wind observations. In § 7 we outline a number of these ideas in a way that should be understandable without a detailed reading of the main text, as well as making more specific predictions that may help to further test and refine the reflection-driven turbulence paradigm. In summary, the model naturally explains the observed decrease in turbulence imbalance with heliocentric radius (Tu & Marsch 1995), as well as its correlation with wind speed if $\chi _{\rm exp}$ is statistically lower in slower streams, as expected from flux-tube expansion models (Chandran 2021). For similar reasons, observations of the classic $\sigma _{c}$-$\sigma _{r}$ circle plot (Bavassano et al. 1998) are reproduced numerically and understandable by appeal to the simple phenomenology and the dominance of long-wavelength structures in the balanced regime. The transition from imbalanced to balanced turbulence also entails a slow shutting off of the turbulent heating, which seems plausibly consistent with observations of the radial and stream dependence of solar-wind heating rates (Totten et al. 1995; Halekas et al. 2023), though more detailed models and observations are needed. Our simulations, as well as previous literature on the subject (Velli et al. 1989; Verdini et al. 2009; Perez & Chandran 2013), reproduce the well-known $1/f$-range spectrum at large scales ($\mathcal {E}^+(k_{\perp })\propto k_{\perp }^{-1}$ in the simulation). Because of the inverse energy transfer, this forms naturally from smaller-scale fluctuations in the initial conditions, migrating to larger scales in the co-moving frame with time. This inverse energy transfer may be observable through its radial dependence or via direct measurements of the turbulent flux, and could have interesting consequences for explaining the dominance of low-frequency fluctuations in observations, even though they should be filtered out by the large Alfvén-speed gradients in the upper chromosphere. Finally, the magnetically dominated Alfvén vortices, which inevitably dominate the solutions at large $a$, seem to bear close resemblance to MFDTs (Tu & Marsch 1991) observed at large heliocentric distances, while the slow (compared with WKB) decay of these structures provides a good fit to observed scalings of fluctuation amplitudes at large heliocentric distances (Zank et al. 1996).

8.2 Uncertainties and future work

Due to both the highly idealized model and the details of the simulation design, our study is beset with a number of significant uncertainties. While we do not believe that these fundamentally invalidate our main results, they are nonetheless important to acknowledge and, hopefully, to rectify in future work.

Setting aside, for a moment, issues with the RMHD EBM itself, the basic phenomenology of § 4.1 (Dmitruk et al. 2002) does not satisfactorily explain some features of the imbalanced-phase turbulence. A priority of future work should be to understand this ‘platonic’ form of reflection-driven turbulence in the expanding box. Of particular difficulty is the relationship between the growth of $\tilde {z}^{-}$, which we observe to be significantly (${\propto }a^{1/2}$) faster than standard treatments (Cranmer & van Ballegooijen 2005; Verdini & Velli 2007; Chandran & Hollweg 2009; van der Holst et al. 2014), and the evolution of the dominant scales of $\tilde {z}^+$ and $\tilde {z}^{-}$ ($\tilde {\lambda }_+$ and $\tilde {\lambda }_{-}$). The growth of $\tilde {\lambda }_+$ and faster-than-expected growth of $\tilde {z}^{-}$ accelerate the transition into the balanced regime, thus decreasing the overall energy decay and heating, so these uncertainties pertain directly to the global energetics of the solar wind, as well as the measured imbalance and residual energy. It may be that some of these discrepancies with the model relate to our initial conditions, and indeed we have found some dependence of the results on the initial conditions (e.g. the infrared spectrum and parallel scales) that remain incompletely understood. Another important goal of future work should be to better explore the dependence on $\varDelta _{\rm box}$, which sets the range of parallel wavelengths available to the system. In our simulations, which fixed $\varDelta _{\rm box}=10$, only the $k_{z}=0$ 2-D mode is linearly expansion dominated (non-Alfvénic), but in reality there should be a continuum of such modes down to the scales where global effects become important ($k_{z}\sim 1/R$). Decreasing $\varDelta _{\rm box}$ is equivalent to increasing the parallel box scale $L_{z}$, and therefore, increases the simulation cost, but this should be explored in future work. An additional priority for future work is to elucidate the physical mechanisms that give rise to the $E^+(k_\perp ) \propto k_\perp ^{-1}$ and $E^-(k_\perp )\propto k_\perp ^{-3/2}$ scalings shown in figure 3(c), which are not explained by existing cascade models for imbalanced MHD turbulence (e.g. Velli et al. 1989; Lithwick et al. 2007; Chandran & Perez 2019).

Moving beyond the uncertainties in interpreting the RMHD EBM results, there exist many uncertainties related to the model itself. Although its simplicity is appealing, RMHD obviously cannot capture any compressive physics or the physics of the large-amplitude spherically polarized fluctuations that are routinely observed in situ (Bale et al. 2019). The latter can be rectified via full expanding-box MHD simulations (Squire et al. 2020; Johnston et al. 2022), but the former arguably cannot, given that the solar wind is a collisionless plasma with compressive fluctuations that may or may not be well described by fluid models (Schekochihin et al. 2009; Verscharen, Chen & Wicks 2017). These issues, as well as our neglect of the Parker spiral, are likely particularly important for our results related to Alfvén vortices, since these structures are inherently compressive (though in total pressure balance). There also exist various subtle issues related to the EBM, motivating future studies with global flux-tube models (van Ballegooijen et al. 2011; Perez & Chandran 2013) that are more focused on super-Alfvénic regions. The EBM should accurately capture dynamics only in the limit where a reflected $\boldsymbol {z}^{-}$ cannot propagate further than one box length, because otherwise this $\boldsymbol {z}^{-}$ could re-encounter the same $\boldsymbol {z}^+$ multiple times (clearly an unphysical effect). This likely limits its applicability to the study of the strong-turbulence regime where $\boldsymbol {z}^{-}$ is anomalously coherent. Another effect that cannot be captured in the EBM due to its fixed parallel size relates to the increased range of long-wavelength, expansion-dominated modes that become available to the system at increasing radius as it transitions into the balanced regime (see § 5.1).

Finally, a key omission, which has been made purely for the sake of simplicity, is the recently discovered ‘helicity barrier’ effect (Meyrand et al. 2021). The helicity barrier suppresses dissipation via electron heating due to finite-Larmor-radius effects in $\beta \lesssim 1$ turbulence, channeling the turbulent flux into ion-cyclotron heating only once the fluctuations can reach a sufficiently small parallel scale (Squire et al. 2022). By suppressing the dissipation of $z^+$, the helicity barrier could significantly change our results in $\beta \lesssim 1$ regions, bringing in direct dependence on the parallel scales. Therefore, our results here can only apply to either the saturated phase, in which the energy flux into ion-gyroradius scales is balanced by ion heating through the cyclotron resonance (Bowen et al. 2024), or to $\beta \gtrsim 1$ regions. Understanding the impact of the helicity barrier on reflection-driven turbulence should be a priority for future work.