2.1. Governing equations for reflection-driven turbulence

The solar wind’s turbulent evolution is governed by the interplay between reflection and nonlinear interactions. To model this in regions beyond the Alfvén critical point, where the flow becomes super-Alfvénic, we adopt the EBM (Grappin et al. Reference Grappin, Velli and Mangeney1993), which approximates the spherical expansion of the solar wind in a Cartesian patch of wind co-moving with the constant radial flow velocity $ U$ . Aligning the $ x$ -axis with the radial direction, the EBM introduces an expansion factor $ a(t) = R(t)/R_0$ , where $ R(t) = R_0 + U t$ is the heliocentric distance for some initial $R_0$ , and modifies spatial gradients as $ \tilde {\boldsymbol{\nabla }} = (\partial _x, a^{-1}\partial _y, a^{-1}\partial _z)$ , viz., $x$ (radial) is unstretched, and ${y},{z}$ (azimuthal/normal) are stretched. Here, $\dot {a} = U/R_0$ is the constant expansion rate due to $U$ . The governing MHD equations in this frame become

(2.1) \begin{align} \frac {\partial \rho }{\partial t} + \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }(\rho \boldsymbol u) &= -2\frac {\dot {a}}{a}\rho , \\[-10pt]\nonumber\end{align}

(2.2) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }}) \boldsymbol{u} &= -\frac {\tilde {\boldsymbol{\nabla }} p}{\rho } + \frac {(\boldsymbol{B} \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }})\boldsymbol{B}}{4\pi \rho } - \frac {\dot {a}}{a} \mathbb{T} \boldsymbol{\cdot }\boldsymbol{u} + \mathcal{D}_{u} , \\[-10pt]\nonumber\end{align}

(2.3) \begin{align} \frac {\partial \boldsymbol{B}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }}) \boldsymbol{B} &= (\boldsymbol{B} \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }}) \boldsymbol{u} - \frac {\dot {a}}{a} \mathbb{L} \boldsymbol{\cdot }\boldsymbol{B} + \mathcal{D}_{B} , \end{align}

where $ \mathbb{T} = \text{diag}(0,1,1)$ and $ \mathbb{L} = \text{diag}(2,1,1)$ encode anisotropic damping from angular momentum conservation and magnetic flux freezing, respectively (here diag denotes a diagonal matrix with the specified values), while $\mathcal{D}_{u}$ and $\mathcal{D}_{B}$ are the dissipative (viscous and resistive) terms that act only at small scales. In the real solar wind they correspond to kinetic damping processes, while in simulations they mainly represent explicit and numerical dissipation near the grid scale.

We assume a locally isothermal equation of state, $p = c_s^2 \rho$ , where $p$ is the thermal pressure, $c_s$ is the sound speed and $\rho$ is the plasma density. By ‘locally isothermal’, we mean that the temperature is spatially uniform throughout the computational domain at any given instant during the simulation. However, this does not imply that the temperature remains constant with heliocentric distance. Rather, we allow the temperature (and thus $c_s$ ) to evolve with $a$ in a manner that reproduces the cooling expected for expanding gas. In particular, the sound speed varies with the expansion as $c_s \propto a^{-2/3}$ , representing the cooling of the solar wind as it expands. The locally isothermal approximation thus captures the background thermal evolution while simplifying the thermodynamics by enforcing spatial temperature uniformity at each time step.

Figure 1.

Schematic of an expanding plasma parcel in the solar wind. (Top) Geometry of the expanding box with axes aligned to radial ( $x$ ), azimuthal ( $y$ ) and normal ( $z$ ) directions. The solar wind initially flows radially outward, while the box expands transversely. The mean magnetic field $\bar {\boldsymbol B}$ (green) is radial near the Sun but rotates into a PS with angle $\varPhi$ as $a(t)$ increases, so that $B_y/B_x \sim a(t)$ . Wave vectors $\boldsymbol{k}$ are shown at angle $\vartheta$ to $\bar {\boldsymbol B}$ , with azimuthal ( $k^{(Y)}$ ) and normal ( $k^{(Z)}$ ) components indicated. (Bottom) Comparison of eddy evolution under (a) purely radial and (b) PS expansion. In the radial case, the perpendicular scale $\ell _\perp$ grows uniformly with expansion. In the PS case, rotation of the mean field changes $\ell _\perp$ creating three-dimensional anisotropy of the eddies.

We take the $y$ -axis to align with the ecliptic so that in the presence of PS the mean magnetic field has both radial ( $x$ ) and azimuthal ( $y$ ) components, ${\bar {\boldsymbol B}} = {B_x} \hat {\boldsymbol x} + B_y\hat {\boldsymbol y}$ . We denote spatial averages over the co-moving expanding domain by $\langle \boldsymbol{\cdot }\rangle$ and fluctuations about a mean by $\delta (\boldsymbol{\cdot })\equiv (\boldsymbol{\cdot })-\langle (\boldsymbol{\cdot })\rangle$ . In particular the mean magnetic field is $\bar {\boldsymbol B}=\langle \boldsymbol{B}\rangle ,$ and we define the corresponding unit vector $\hat {\boldsymbol b}=\bar {\boldsymbol B}/\lvert {B}\rvert .$ From $\hat {\boldsymbol b}$ we form an orthonormal triad $(\hat e_\parallel ,\hat e_T,\hat e_N$ , see figure 1) by taking

(2.4) \begin{align} \hat e_\parallel &= \hat {\boldsymbol b}, \quad \hat e_T = \frac {\hat {\boldsymbol z}\times \hat {\boldsymbol b}}{\lvert \hat {\boldsymbol z}\times \hat {\boldsymbol b}\rvert }, \quad \hat e_N = \hat e_\parallel \times \hat e_T=\hat {\boldsymbol z}. \end{align}

The magnetic-field components in this triad are

(2.5) \begin{equation} B_\parallel =\boldsymbol{B}\boldsymbol{\cdot }\hat e_\parallel ,\qquad B_T=\boldsymbol{B}\boldsymbol{\cdot }\hat e_T,\qquad B_N=\boldsymbol{B}\boldsymbol{\cdot }\hat e_N. \end{equation}

Due to expansion, the plasma mean quantities evolve with time. The magnetic-field components follow $B_x = B_{x0}/a^2$ , $B_y = B_{y0}/a$ and $B_z = B_{z0}/a$ , and the density as $\rho = \rho _0/a^2$ , where, the subscript 0 is to refer to the values at $t = 0 (a = 1)$ . The Alfvén velocity becomes

(2.6) \begin{equation} \boldsymbol{v}_{\textrm {A}} \equiv \frac {{\bar {\boldsymbol B}}}{\sqrt {4\pi \rho }} = \frac {v_{\textrm {A0},\textit {x}}}{a} \hat {\boldsymbol x} + {v}_{\textrm {A0},\textit {y}} \hat {\boldsymbol y}, \quad {v}_{\textrm {A0},\textit {x}} = \frac {B_{x0}}{\sqrt {4\pi \rho _0}}, \quad {v}_{\textrm {A0},\textit {y}} = \frac {B_{y0}}{\sqrt {4\pi \rho _0}}, \end{equation}

such that the radial component ${v}_{\textrm {A},\textit {x}} \propto 1/a$ decreases with expansion, while ${v}_{\textrm {A},\textit {y}}$ remains constant. Based on these introduced scalings, the Parker angle $\varPhi$ , which quantifies the deviation of the background magnetic field from the radial direction in the $x-y$ -plane, is

(2.7) \begin{equation} \tan \varPhi = \frac {\langle {B}_y\rangle }{\langle {B}_x\rangle } \; =\; a \frac {{B}_{y0}}{{B}_{x0}}\propto a, \end{equation}

indicating that, if the initial background magnetic field has a non-zero transverse component, it will increasingly deviate from the radial direction as the system undergoes expansion. This is the manifestation of PS field in the local domain.

To develop our phenomenology, we assume locally constant density $ \rho$ and incompressible flow ( $ \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{u} = 0$ ). This allows us to simplify the equations and focus on the predominantly Alfvénic fluctuations. To study the evolution of these Alfvén waves – particularly the important phenomenon of reflection caused by gradients in the mean Alfvén speed – it is standard to reformulate the MHD equations in terms of the Elsässer variables: $ \boldsymbol{z}^\pm = \boldsymbol{u} \mp \boldsymbol{B}/\sqrt {4\pi \rho } = \boldsymbol{u} \mp \boldsymbol{b} \mp \boldsymbol{v_{\mathrm{A}}}$ (where $\boldsymbol{b} = \delta \boldsymbol{B}/\sqrt {4\pi \rho }$ ). The Elsässer formalism is useful because $\boldsymbol z^+$ and $\boldsymbol z^-$ signs distinguish fluctuations that propagate parallel and anti-parallel to the background field, respectively, clarifying the interplay of counter-propagating Alfvénic fluctuations. Combining the momentum equation (2.2) and the induction equation (2.3) yields

(2.8) \begin{align} \frac {\partial \boldsymbol{z}^\pm }{\partial t} \pm \big( \boldsymbol{v}_A \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \big) \boldsymbol{z}^\pm + \big( \boldsymbol{z}^\mp \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \big) \boldsymbol{z}^\pm + \tilde {\boldsymbol{\nabla }} {p \prime }^\pm = -\frac {\dot {a}}{2a} \mathbb{T} \boldsymbol{\cdot }\left ( \boldsymbol{z}^+ + \boldsymbol{z}^- \right ) \nonumber \\ - \frac {\dot {a}}{2a} \left ( z_x^\pm - z_x^\mp \right ) \hat {\boldsymbol x} + \mathcal{D}^\pm , \end{align}

where $\mathcal{D}^\pm$ represents the dissipative effects due to viscosity/resistivity that act on $\boldsymbol z^\pm$ , and $\tilde{\boldsymbol{\nabla}}p'^\pm$ ensures $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{z}^\pm = 0$ . We take the mean magnetic field $\bar {\boldsymbol B}$ to point in the positive radial direction, so that $\boldsymbol z^+$ perturbations propagate outward in the absence of reflection.

2.2. Wave-action Elsässer variables

In an expanding medium like the solar wind, Alfvén wave amplitudes can naturally decay due to geometric Wentzel–Kramers–Brillouin (WKB) effects as well as true dissipation (Heinemann & Olbert Reference Heinemann and Olbert1980). To remove this trivial evolution, which manifests in (2.8) via $\boldsymbol z^\pm$ terms on right-hand side, we rescale the Elsässer variables as (Chandran & Hollweg Reference Chandran and Hollweg2009; Meyrand et al. Reference Meyrand, Squire, Mallet and Chandran2025)

(2.9) \begin{equation} \tilde {\boldsymbol z}^\pm \;=\;a^{1/2}\,\boldsymbol z^\pm \;\propto \;\frac {\boldsymbol z^\pm }{\sqrt {\omega _{\textrm {A}}}}, \end{equation}

where $\omega _{\textrm {A}} = \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{v}_{\textrm {A}} = k_{x,0} (v_{\textrm {A0},\textit {x}}/a) + (k_{y,0}/a) v_{\textrm {A0},\textit {y}} \propto 1/a$ is the local Alfvén frequency for wavenumber $\boldsymbol k$ . This transformation isolates the physical processes – nonlinear cascade and reflection – that genuinely redistribute and dissipate wave energy. The nearly conserved quantity in this framework is the WKB wave-action density $|\boldsymbol z^\pm |^2/\omega _{\textrm {A}}$ , which remains constant without reflection or nonlinear coupling. Rewriting (2.8) in terms of $\tilde {\boldsymbol z}^\pm$ cancels all expansion–induced decay, leaving

(2.10) \begin{equation} \dot {a}\frac {\partial \tilde {\boldsymbol{z}}^\pm }{\partial a} \pm \big( \boldsymbol{v}_A \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \big) \tilde {\boldsymbol{z}}^\pm + {a^{-1/2}} \big( \tilde {\boldsymbol{z}}^\mp \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \big) \tilde {\boldsymbol{z}}^\pm + \tilde {\boldsymbol{\nabla }} \tilde {p}' = -\frac {\dot {a}}{2a} \tilde {\mathbb{T}}\boldsymbol{\cdot }\tilde {\boldsymbol{z}}^\mp + \mathcal{\tilde {D}}^\pm , \end{equation}

where $\tilde {\mathbb{T}}= \textrm {diag} (-1,1,1)$ . Equation (2.10) shows explicitly how reflection seeds $\boldsymbol z^-$ from $\boldsymbol z^+$ fluctuations, while nonlinear interactions, which weaken in time due to the $a^{-1/2}$ factor, can sustain a cascade. The final term on the right-hand side of (2.10), $-(\dot {a}/2a)\tilde {\mathbb{T}}\boldsymbol{\cdot }\tilde {\boldsymbol{z}}^\mp$ , represents an anisotropic reflection effect associated with the expansion. In the radial field case with perpendicular fluctuations, $\tilde z^\pm _x$ is negligible and the dominant reflection occurs through $\tilde z_{y,z}^- \propto -\tilde z_{y,z}^+$ , thereby naturally causing $\langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle \lt 0$ ; where $\langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle /2 = (\langle \boldsymbol u^2\rangle -\langle \boldsymbol b^2\rangle )/2$ is the wave-action residual energy $\tilde E_{r}$ . However, with either large-amplitude spherically polarised fluctuations or an oblique background field, there will be a significant $z_x^\pm$ component, so that we also have contribution from $\tilde z_x^- \propto \tilde z_x^+$ . This anisotropic coupling could thus moderate the growth of the cross-correlation $\langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle$ , resulting in systematically smaller values of $\lvert \langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle \rvert$ in the PS configuration compared with the radial case. Physically, this implies that strong negative correlations between counter-propagating Elsässer fields are expected in the radial case, whereas in the PS geometry the anisotropic reflection term reduces the magnitude of $\langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle$ .

2.3. Energy dynamics in expanding solar-wind turbulence

In MHD turbulence, energy is distributed between kinetic and magnetic fluctuations, with nonlinear interactions transferring energy across scales. The wave-action Elsässer energy densities are defined by $\tilde {E}^\pm \equiv \langle |\tilde {\boldsymbol{z}}^\pm |^2 \rangle / 4$ . Unlike homogeneous MHD these energies are not ideally conserved under solar-wind expansion, because reflection acts as either a source or sink depending on the correlation between counter-propagating modes. Multiplying the evolution (2.10) by $\tilde {\boldsymbol{z}}^\pm$ , then averaging gives

(2.11) \begin{equation} \dot {a} \frac {\partial }{\partial a} {\tilde {E}^\pm } = -\frac {\dot {a}}{4a} { \langle \tilde {\boldsymbol{z}}^+ \boldsymbol{\cdot }\mathbb{T} \boldsymbol{\cdot }\tilde {\boldsymbol{z}}^- \rangle } + \frac {\dot {a}}{4a} \langle \tilde {z}_x^+ \tilde {z}_x^- \rangle + \mathcal{\tilde {D}}^\pm , \end{equation}

because we expect predominantly Alfvénic fluctuations, which are perpendicular to $\bar {\boldsymbol B}$ , in the case of the purely radial field the term $\langle \tilde {z}_x^+ \tilde {z}_x^- \rangle$ is likely negligible, so the energy evolution simplifies to $\dot {a} \, \partial \tilde {E}/\partial a = -(\dot {a}/2a) \tilde E_{r}$ (Meyrand et al. Reference Meyrand, Squire, Mallet and Chandran2025). The reflection thus acts as a source or sink of wave-action energy depending on the sign of the residual energy $\tilde E_{r}$ (note that any change in fluctuation energy is balanced by the background-flow energy budget; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015). In the absence of dissipation and external forcing, $\tilde {E}^+ - \tilde {E}^-$ is conserved by the combined dynamics of wave reflection and advection, because the corresponding source terms are identical in the $+$ and $-$ equations and cancel upon subtraction; making this difference the true wave-action invariant (Heinemann & Olbert Reference Heinemann and Olbert1980).

In the PS configuration, the $+ \langle \tilde {z}_x^+ \tilde {z}_x^- \rangle$ term could become significant even when $\boldsymbol z^\pm$ is predominantly perpendicular to $\boldsymbol B$ , and we likewise see from (2.10) that it is naturally driven with the correlation required to act as a source of $\tilde {E}^\pm$ .

2.4. Reflection-driven turbulence phenomenology

The RDT phenomenology involves the consideration of different balance of key terms. It is helpful to sketch (2.10) as

(2.12) \begin{equation} \dot {a}\frac {\partial \tilde {\boldsymbol{z}}^\pm }{\partial a} + \mathcal{L}^\pm \tilde {\boldsymbol{z}}^\pm + \mathcal{N}^\pm \tilde {\boldsymbol{z}}^\pm = - \mathcal{R}^\pm \tilde {\boldsymbol{z}}^\mp , \end{equation}

where $\mathcal{L}^\pm = \boldsymbol{v}_A \boldsymbol{\cdot }\tilde{\boldsymbol{\nabla}}$ represents Alfvénic propagation, $\mathcal{N}^\pm = {\boldsymbol{z}}^\mp \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} = a^{-1/2} \tilde {\boldsymbol z}^\mp \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \sim {z}^\mp /\ell _\perp$ represents nonlinear terms that will damp out ${z}^\pm$ (here, $\ell _\perp$ denotes the perpendicular correlation length scale of fluctuations and ${z}^+$ denotes the characteristic outer-scale amplitude of outward-propagating fluctuations at scale $\ell _\perp$ ) and reflection is $\mathcal{R}^\pm = -\dot {a} / 2 a \mathbb{\tilde T}$ . The characteristic time scales associated with the terms in (2.12) can differ in their scaling. The Alfvén propagation time along the mean field is $\tau _{\textrm {A}} = \ell _\parallel /v_{\textrm {A}}$ ( $\ell _\parallel$ is the parallel correlation length), while the nonlinear interaction time across the field is $\tau _{\textrm {nl}} = \ell _\perp /{z}^\mp$ . Both $\tau _{\textrm {A}}$ and $\tau _{\textrm {nl}}$ decrease at smaller scales, since they depend directly on the characteristic sizes of turbulent structures along and across the mean magnetic field, respectively. In contrast, the expansion time $\tau _{\textrm {exp}}= a/\dot {a}$ is independent of scale. Following Meyrand et al. (Reference Meyrand, Squire, Mallet and Chandran2025), we define the dimensionless ratios

(2.13) \begin{equation} \chi _{\textrm {exp}} = \frac {\mathcal{N} }{\mathcal{R}} = \frac {\tau _{\textrm {exp}}}{\tau _{\textrm {nl}}}, \quad \chi _{\textrm {A}} = \frac {\mathcal{N} }{\mathcal{L}} = \frac {\tau _{\textrm {A}}}{\tau _{\textrm {nl}}}, \end{equation}

which quantify the relative importance of expansion and Alfvén propagation compared with nonlinear interactions (Goldreich & Sridhar Reference Goldreich and Sridhar1995). These ratios will play a key role in controlling the dynamics of the turbulent cascade in an expanding solar wind.

The standard phenomenology (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Verdini & Velli Reference Verdini and Velli2007; Chandran & Hollweg Reference Chandran and Hollweg2009) balances these terms to model turbulence and heating, making two key balancing assumptions. First, for $\tilde {z}^-$ (backward waves) propagation and time variation are assumed negligible, the reflection source ( $\mathcal{R}^- \tilde {z}^+$ ) balances the nonlinear sink ( $\mathcal{N}^- \tilde {z}^-$ )

(2.14) \begin{equation} \mathcal{N}^- \tilde {z}^- \sim \mathcal{R}^- \tilde {z}^+\implies \tilde {z}^- \sim \frac {\mathcal{R}^- \tilde {z}^+}{\mathcal{N}^-} \sim \frac {\dot {a}}{2a}\frac {\tilde {z}^+}{({ z}^+/\ell _\perp )} = a^{1/2}\left (\frac {\dot {a}}{2a}\right ) \ell _\perp , \end{equation}

and using $(\dot a/a)/(z^+ /\ell _\perp )=\chi _{\textrm {exp}}^{-1}$ gives

(2.15) \begin{equation} \tilde {z}^- \simeq \frac {\tilde {z}^+}{\chi _{\textrm {exp}}}\; \Rightarrow \; \frac {{z}^-}{{z}^+} \simeq \chi _{\textrm {exp}}^{-1} . \end{equation}

This means that the $\tilde {z}^-$ amplitude is small – i.e. turbulence is highly imbalanced – when reflection is weak or nonlinearity is strong $\chi _{\textrm {exp}} \gg 1$ . The second assumption concerns the decay of $\tilde {\boldsymbol z}^+$ . In the $\tilde z^+$ equation we neglect reflection because $\tilde z^-\ll \tilde z^+$ . The linear propagation terms ( $\mathcal L^\pm \tilde {\boldsymbol z}^\pm$ ) does not directly lead to decay – it merely transports fluctuations along the mean field – so $\partial \tilde {\boldsymbol z}^+/\partial t$ is dominated by damping at the nonlinear rate ( $\mathcal N^+$ )

(2.16) \begin{equation} \dot {a}\frac {\partial \tilde {{z}}^+}{\partial a} \simeq - \mathcal{N}^+ \tilde {{z}}^+ \implies \dot {a}\frac {\partial }{\partial a}\tilde {{z}}^+ \simeq - \frac {{ z}^-}{\ell _\perp }\tilde {{z}}^+ \simeq -\frac {\dot {a}}{2 a} \tilde {z}^+ \implies \tilde {z}^+ \propto a^{-1/2}. \end{equation}

This corresponds to $z^+\propto a^{-1}$ or (in a radial field) $z^+/v_{\textrm {A}} \sim \text{const}$ . The heating rate $\mathcal{Q}$ , which is the rate at which energy in the outward-propagating mode $\tilde {z}^+$ is damped, can be calculated using (2.11) by multiplying with $\rho V$ to get actual energy (see, Perez et al. Reference Perez, Chandran, Klein and Martinović2021)Footnote ¹ , where $ V \propto a^2$ is the volume of the domain. This gives

(2.17) \begin{equation} \mathcal{Q}= \frac {\dot {a}}{a} \frac {\partial }{\partial a} {(\rho \tilde {E}^+ a^2)} = \rho _0\frac {\dot {a}}{a} \frac {\partial }{\partial a} \frac {\langle |\tilde {\boldsymbol z}^+|\rangle ^2}{4} \approx \rho _0 \frac {\dot {a} \tilde {z}^+}{2a} \frac {\partial \tilde {z}^+}{\partial a} \approx \rho _0\frac {\dot {a}}{4 a^2} (\tilde {z}^+)^2, \end{equation}

where we neglected the residual energy and use $\langle |\tilde {\boldsymbol z}^+|\rangle ^2 \approx (\tilde {z}^+)^2$ (mean square amplitude). The phenomenology in § 2.3 predicts $\tilde {E}\propto a^{-1}$ , which implies a heating rate per co-moving volume $\mathcal{Q}\propto a^{-3}$ . Converting to physical units (physical volume $\propto a^2$ ) gives a heating rate per physical volume $\mathcal{\tilde {Q}}\propto a^{-5}$ . The heating rate $\mathcal{Q}$ represents the actual dissipation of turbulent energy into thermal energy through viscous and resistive processes, not the reversible redistribution of energy between outward- and inward-propagating modes. Reflection converts $\boldsymbol z^+$ into $\boldsymbol z^-$ and vice versa, and the nonlinear turbulent cascade transfers energy from large scales to small scales, but only when the cascade reaches dissipative scales is energy irreversibly removed and counted in $\mathcal{Q}$ . In our simulations, explicit viscosity and resistivity are not included; instead, numerical dissipation at the grid scale in Athena++ serves as an effective dissipation mechanism that captures the cascade. The rate $\mathcal{Q}$ therefore measures how rapidly the outward energy $\tilde {E}^+$ is depleted by this dissipative process, which is fundamentally distinct from the reversible energy exchange between $\boldsymbol z^+$ and $\boldsymbol z^-$ driven by reflection. We note that while kinetic effects likely dominate dissipation in the actual solar-wind plasma, our MHD framework with numerical dissipation at the grid scale captures the essential phenomenology of energy removal at small scales.

Interestingly, this heating rate does not depend on the nonlinear rate or scale of the turbulence, because a smaller transverse scale $z^+$ damps $z^-$ more rapidly, a smaller $\ell _\perp$ implies that $z^-$ has smaller amplitude for the same reflection rate $\mathcal{R}$ . This lower amplitude then cancels the increased nonlinear rate coming from the smaller scale of $z^-$ , thus leading to the same overall efficiency at damping $z^+$ .

The RDT framework rests on following assumptions that demand careful scrutiny. First, the nonlinear damping rate provided by the turbulence is $ {z}^\mp /\ell _\perp$ , implying strong turbulence ( $\chi _{\textrm {A}} \sim 1$ ) for validity (but see related weak turbulence version in Chandran & Perez (Reference Chandran and Perez2019)), which leads to the same scaling. Second, the phenomenology assumes that a single dominant scale $\ell _\perp$ governs both $\tilde {z}^+$ and $\tilde {z}^-$ ; this $\ell _\perp$ cancels in the heating rate $\mathcal{Q}$ , but in principle could be different for $\tilde {z}^+$ and $\tilde {z}^-$ (indeed, this was seen by Meyrand et al. (Reference Meyrand, Squire, Mallet and Chandran2025)). Third, anomalous coherence, with $\tilde {z}^-$ carried along by $ \tilde {z}^+$ , is required to preserve the $\tilde {z}^-$ balance and to justify neglecting the ( $ v_{\textrm {A}} \boldsymbol{\cdot }\boldsymbol{\nabla}$ ) term in our first assumption. Fourth, $\tilde {z}^-$ must adjust to reflection-nonlinear equilibrium faster than background changes; this requirement reduces to $\chi _{\textrm {exp}} \gtrsim 1$ . Fifth, reflection in the $\tilde {z}^+$ equation must be negligible ( $\mathcal{R}^+ \tilde {z}^- \ll \mathcal{N}^+ \tilde {z}^+$ ), which likewise reduces to the same $\chi _{\textrm {exp}} \gtrsim 1$ condition.

Thus the phenomenology is self-consistent only for sufficiently large amplitude $\tilde z^+$ such that $\chi _{\textrm {exp}}\gg 1$ , which implies strong imbalance $\tilde z^-\ll \tilde z^+$ . We therefore expect the system to transition to the balanced regime once the outward Elsässer energy decays to the point that $\chi _{\textrm {exp}}\sim 1$ . Conversely, when $\chi _{\textrm {exp}} \lt 1$ , the ordering breaks down, reflection becomes dominant and the dynamics from (2.17) likely reverts toward reflection-dominated (essentially linear) behaviour. In this regime with a radial $\bar {\boldsymbol B}$ , the system evolves into a magnetically dominated, balanced state in which $\tilde {\boldsymbol z}^+ \simeq -\tilde {\boldsymbol z}^-$ , and $\tilde {E}$ starts growing in time (see (2.11)). Nonlinear dissipation is strongly suppressed, so the turbulent heating effectively shuts off (Meyrand et al. Reference Meyrand, Squire, Mallet and Chandran2025).

2.4.1. The role of Parker spiral

The key new contribution of this work is the extension of RDT phenomenology to the PS geometry. We argue that the underlying dynamics remains similar to that in a radial field, but that the spiral geometry introduces a key geometric modification affecting the evolution of eddy shapes. In a purely radial field, as the solar wind expands, eddies are stretched by expansion preferentially transverse to the magnetic-field direction, becoming increasingly pancake shaped: $\ell _\perp$ increases while $k_\perp$ decreases. This stretching rapidly reduces the nonlinear turnover rate, accelerating the approach to $\chi _{\textrm {exp}}= 1$ . However, as the PS twists the mean magnetic field away from purely radial, it effectively rotates the mean field relative to how eddy structures are stretched by expansion. This rotation causes the field to ‘cut across’ different parts of the eddies, so the eddies no longer remain perfect pancakes with respect to the local mean field and develop three-dimensional, anisotropic outer-scale shapes. A cartoon of this physics is shown in figure 1. To capture this anisotropic deformation, we use the local orthogonal basis $(\hat {e}_\parallel , \hat {e}_T, \hat {e}_T)$ to define two transverse outer scales, $\ell _{\perp ,\textrm {T}}$ and $\ell _{\perp ,\textrm {N}}$ , characterising the turbulent eddy scales along $\hat {e}_T$ and $\hat {e}_N$ . The nonlinear interaction time is then expected to be determined by the smaller of these two transverse scales

(2.18) \begin{equation} \tau ^\pm _{\textrm {nl}} \sim \frac {\min (\ell _{\perp ,\textrm {T}},\; \ell _{\perp ,\textrm {N}})}{{z}^\mp }. \end{equation}

As $\ell _\perp (a)$ starts flattening or even decreases with $a$ , this decreases the outer scale’s nonlinear time scale $\tau _{\textrm {nl}}$ . Consequently, this helps maintain $\chi _{\textrm {exp}}$ large enough to prolong the imbalanced cascade. In this way, the competition between rotation and expansion means that the PS geometry prevents the eddies from becoming indefinitely flattened into pancakes, sustaining turbulence and heating farther out in the solar wind.

To quantify that intuition, we will calculate – under the simplifying assumption that outer scales evolve linearly with the background flow – how the wave obliquity, projected wave vector angles, perpendicular wavenumber $k_\perp$ , nonlinear time $\tau _{\textrm {nl}}$ and thus $\chi _{\textrm {exp}}$ evolve with expansion. For the evolution of $z^+$ , we adopt the same phenomenology discussed above, with $\tilde {z}^+\propto a^{-1/2}$ . This is motivated by the fact that the reflection, although differing in sign in the $x$ -direction, has the same magnitude in all directions. Related effects were previously explored in the context of switchbacks by Squire et al. (Reference Squire, Johnston, Mallet and Meyrand2022). It is helpful to define the following angles: the wave obliquity parameter $\vartheta \approx \arccos {(\hat {\boldsymbol{k}} \boldsymbol{\cdot }\bar {\boldsymbol B})}$ quantifies the angle between the wave vector direction $\hat {\boldsymbol{k}}$ and the mean magnetic field $\bar {\boldsymbol B}$ , which determines $\ell _\perp$ ; $\theta _{p0}$ measures the initial angle between the wave vector and the radial direction; and $\varphi$ denotes the angle of the wave vector projected onto the tangential–normal ( $y$ – $z$ ) plane, such that $\varphi = 0$ corresponds to $\boldsymbol{k}$ aligned with the tangential direction and $\varphi = \pi /2$ to $\boldsymbol{k}$ aligned with the normal direction (figure 1).

Figure 2.

Evolution of key parameters for waves under solar-wind expansion, plotted versus expansion factor $a$ . Thick black: radial field ( $\varPhi _0=0^\circ$ ). Coloured: PS cases ( $\varPhi _0=2^\circ{-}20^\circ$ , darker colours correspond to smaller $\varPhi _0$ ). (a) Wave amplitude $z^+/v_{\textrm {A}}$ stays roughly constant for the purely radial case, and in the PS case remains constant initially but then decays as $\propto a^{-1}$ once the azimuthal component becomes significant and $v_{\textrm {A}}\propto a^0$ . (b) Obliquity $\sin \vartheta (a)$ : the angle between $\boldsymbol k$ and $\bar {\boldsymbol B}$ initially decreases and then increases again as the mean field rotates azimuthally, producing a clear inflection with a PS. (c) Expansion cascade parameter $\chi _{\textrm {exp}} = (k_\perp \,z^+_\perp )/(\dot a/a)$ for the out-of-plane case ( $\varphi = \pi /2$ ): turbulence is sustained while $\chi _{\textrm {exp}} \gtrsim 1$ ; $\chi _{\textrm {exp}}$ initially decays as $\propto a^{-1}$ for both radial and PS cases, but for PS it starts increasing and eventually flattens as $a^{0}$ once the azimuthal component becomes dominant. (d) $\chi _{\textrm {exp}}$ for in-plane wave vectors: shown for $\varphi =\pi$ (solid) and $\varphi =0$ (dashed); in both cases, the wave vector lies in the $x$ – $y$ plane. For $\varphi =0$ , the wave passes through purely parallel propagation.

As shown in figure 2, for the purely radial magnetic-field case ( $\varPhi _0=0$ ), the system expands such that the eddy scale grows ( $\ell _\perp \propto a$ ), the fluctuation amplitude decays ( $z^+ \propto a^{-1}$ ) and $v_{\textrm {A}}\propto a^{-1}$ ; the ratio $ z^+/v_{\textrm {A}}$ is thus constant with expansion. In the PS case the growing azimuthal component causes $v_{\textrm {A}}$ to decrease more slowly (and eventually become nearly constant $v_{\textrm {A}} \propto a^{0}$ ), so $z^+/v_{\textrm {A}}$ first remains constant and then decreases approximately as $a^{-1}$ once the spiral angle becomes large. The nonlinear time scale increases because $k_\perp \propto a^{-1}$ , giving $k_\perp z^+ \propto a^{-2}$ ( $\tau _{\textrm {nl}} \propto a^2$ ). This implies $\chi _{\textrm {exp}}$ also decays as $\chi _{\textrm {exp}} \propto a^{-1}$ , rapidly approaching $\chi _{\textrm { exp}}\leqslant 1$ , where the phenomenology breaks.

In the PS case, $\sin \vartheta$ initially decreases as in the radial case, but as $a$ approaches $a_{\min } = \sqrt {\cot \theta _{p0} \,\cot \varPhi _0}$ , its evolution deviates: $\sin \vartheta$ reaches a minimum and then increases again as $\overline {\boldsymbol B}$ rotates toward the perpendicular direction (see figure 2(b) and also Squire et al. Reference Squire, Johnston, Mallet and Meyrand2022). While the detailed dynamics is complex and nonlinear, we assume the overall RDT framework still applies such that $z^+\propto a^{-1}$ , but the evolving geometry changes the scaling of $k_\perp$ . Initially $k_\perp$ decays with expansion as in the radial case, but as the mean field rotates relative to the eddies, making wave vectors oblique and driving $k_\perp$ back up, thereby increasing the nonlinear turnover rate. As a result, instead of continuously decaying, $\chi _{\textrm {exp}}$ levels off staying roughly constant ( $\propto a^0$ ) or even growing at intermediate times, delaying indefinitely the turbulence shutoff that occurs at $\chi _{\textrm {exp}} \lt 1$ and maintaining the high imbalance (see figure 2).

This revival is complicated by wave orientation in the $\hat {e}_T$ – $\hat {e}_N$ plane. When wave vectors lie perpendicular to the spiral plane $(\varphi = \pi /2, \; p^{(Z)})$ , $\chi _{\textrm {exp}}$ exhibits a robust rebound. For in-plane orientations $(p^{(Y)})$ , however, waves can pass through a purely parallel propagation phase $(\vartheta \to 0)$ suggesting a temporary collapse of nonlinear interactions before partial recovery. However, this is likely not significant: it occurs only when the wave vectors lie perfectly in the plane of the PS ( $\varphi =0$ ), which is a very special case. As discussed above, the PS geometry squeezes eddies into three-dimensional anisotropic structures with the shortest perpendicular scale likely determining the nonlinear turnover time. Thus, even $\ell _\perp$ becomes large in one direction the cascade can presumably continue. Note that an important caveat about the assumption $\ell _y=\ell _z \propto a$ may not hold strictly, since MHD nonlinear interactions tend to elongate structures preferentially along the local field direction $\hat {e}_\parallel$ , thereby presumably changing the perpendicular growth for the linear evolution discussed above.

We now test the above ideas with a controlled numerical experiment suite. Using compressible expanding-box MHD simulations, we evolve initially outward $z^+$ fields and measure the properties introduced above to validate the phenomenology and quantify its parameter dependence. The following sections present these numerical results, compare them directly with the theoretical expectations, and highlight where the simulations confirm, refine or challenge our simple model.

3.1. Expanding-box model numerical

To solve the EBM (2.1)–(2.3), we use the finite-volume astrophysical code Athena++ (Stone et al. Reference Stone, Tomida, White and Felker2020). We employ the Harten–Lax–van Leer discontinuities Riemann solver (Mignone Reference Mignone2007), modified to include expansion effects. To simplify the expansion terms and enhance numerical stability, we transform variables as follows (Johnston et al. Reference Johnston, Squire, Mallet and Meyrand2022):

(3.1) \begin{equation} \rho = \lambda ^{-1} \rho ', \quad \boldsymbol{u} = \varLambda \boldsymbol{u}', \quad \boldsymbol{B} = \lambda ^{-1} \varLambda \boldsymbol{B}', \quad \boldsymbol{\nabla }' = \varLambda \tilde {\boldsymbol{\nabla }}, \end{equation}

where $\varLambda =\mathrm{diag}(1,a,a)$ and $\lambda =a^2$ . Following this variable change and expressing the EBM MHD equations in conservative form, we obtain

(3.2) \begin{equation} \frac {\partial \rho '}{\partial {t}} + {\boldsymbol{\nabla }}' \boldsymbol{\cdot }(\rho ' \boldsymbol{u}') = 0, \end{equation}

(3.3) \begin{equation} \frac {\partial \left (\rho '\boldsymbol{u}'\right )}{\partial t} + {\boldsymbol{\nabla }}' \boldsymbol{\cdot }\left (\rho ' \boldsymbol{u}' \boldsymbol{u}' + \left (p+\frac {B^2}{2}\right ) \varLambda ^{-2} - \frac {\boldsymbol{B}' \boldsymbol{B}'}{\lambda }\right ) = -\left ( \frac {d \ln {\lambda }}{dt}\right ) \mathbb{T}'\boldsymbol{\cdot }\rho '\boldsymbol{u'},\end{equation}

(3.4) \begin{equation} \frac {\partial \boldsymbol{B}'}{\partial {t}}-\boldsymbol{\nabla }'\times \left (\boldsymbol{u}'\times \boldsymbol{B}'\right )= 0. \end{equation}

Equations (3.2)–(3.4) are quite similar to the ideal MHD equations, since the continuity and induction equations have no expansion or source terms. Instead, the momentum equation and the variable defined in (3.1), now incorporates all the effects of expansion.

3.2. Simulations set-up and initial conditions

Numerical details of all simulations are given in table 1. All simulations employ periodic boundary conditions in $x,y,z$ , as is standard for local EBM studies. This is justified because the computational patch remains small compared with the heliocentric distance, allowing the EBM mapping of local spherical expansion onto a Cartesian, periodic domain. Periodicity reduces computational cost and therefore enables higher numerical resolution to capture the turbulent cascade across a broad range of scales. We stress the limitations: periodic domains omit open-boundary effects (global inflow/outflow) and can box limit the longest parallel modes (affecting $\ell _\parallel$ and near-two-dimensional structures). To mitigate these effects we use an initially anisotropic (elongated) domain so transverse stretching is naturally represented as $L_y,L_z\propto a$ . Within these bounds the model effectively captures the local, expansion-driven dynamics of RDT beyond the Alfvén critical point.

Table 1.

Simulation parameters for the expanding-box runs. High-resolution (HR) runs use $1200^3$ grid points; moderate-resolution (MR) runs use $720^3$ . The parameters explored are; $\varPhi _0$ the initial PS angle, $L_{x0}/L_{\perp 0}$ the initial box aspect ratio (where $L_{\perp 0}$ denotes both $L_{y0}$ and $L_{z0}$ ), $\dot a$ the expansion rate, $A_0$ the initial fluctuation amplitude and $\beta _0$ the plasma beta. Also, ${k}_{\textrm {peak}}$ sets the centre of the Gaussian peak in $k$ -space, $\chi _{\textrm {exp,0}}$ the initial expansion-to-nonlinearity ratio and $k_{\textrm {w}}$ sets the width of the Gaussian peak.

We set the initial box dimensions to $L_{x0}=1$ and $L_{y0}=L_{z0}=0.2$ . We use the expansion factor $a$ to parameterise the evolution of the system, where $a$ is directly related to heliocentric distance via $a=R/R_0$ , with $R_0$ being the radial distance at which the simulation is initialised ( $a=1$ ). The choice of $R_0$ is arbitrary in the sense that it sets the reference scale and simply rescales other quantities (e.g. $v_{\textrm {A}}$ ) without changing the dimensionless physics. To aid interpretation and correspondence with observational regimes, we can map $a$ to actual solar distances using the evolving PS angle $\varPhi (a)$ . We assume a PS angle of $45^\circ$ occurs at 1 AU (Borovsky Reference Borovsky2010). In our simulations with initial PS angle $\varPhi _0=2^\circ$ , the value $a\, \approx \,28.64$ corresponds to 1 AU, while for $\varPhi _0=5^\circ$ , we have $a\approx 11.5$ at 1 AU. Importantly, the mapping between the heliocentric distance and PS is not unique: it depends on the PS calibration (e.g. the reference angle at 1 AU), the latitudes and on how one chooses the simulation inner boundary. At high latitudes, the PS is much weaker, so off-ecliptic regions remain nearly radial even out to several AU. Note that scanning over $\varPhi _0$ at fixed initial conditions is best interpreted as starting from the same initial $\chi _{\textrm {exp,0}}$ (discussed below) at different heliocentric distances.

As the system expands, the box elongates in the transverse directions; by $a=50$ , the domain becomes strongly oblate, resembling a pancake-like structure. At $a=1$ , we initialise the background magnetic field $\bar {\boldsymbol B}$ in the $x-y$ plane with $|{\bar {\boldsymbol B}}|=1$ and $B_{y,0} \lt 0$ , defining the initial PS angle $\varPhi _0$ .

Waves exhibiting nearly constant magnetic-field strength and strong Alfvénic correlation $(\delta \boldsymbol B_\perp \propto \delta \boldsymbol u_\perp )$ are observed in the solar wind. These waves are called spherically polarised because the constant $|\boldsymbol B|$ is preserved in the fluctuations. This property can be quantified through

(3.5) \begin{equation} C_{B^2} \equiv \frac {\delta \left (| B|^2\right )}{\left (\delta \boldsymbol{B}\right )^2}={\frac { \big({{ \big\langle \big(B^2 - \langle B^2\rangle \big)^2\big\rangle }}\big)^{1/2}}{\left \langle |\boldsymbol{B} - \langle \boldsymbol{B}\rangle |^2\right \rangle }}, \end{equation}

which measures how the components of $\boldsymbol B$ are correlated to keep $|\boldsymbol B|$ constant; $C_B^2=0$ for perfectly spherically polarised fluctuations. The initial condition is chosen to mimic a strongly $z^+$ -dominated, nearly spherically polarised Alfvénic wind propagating outwards from near the Alfvén point: $\boldsymbol z^+(t=0)$ carries the fluctuation energy while $\boldsymbol z^-(t=0)=0$ . Radial expansion and large-scale inhomogeneity produce partial reflection (creation of $\boldsymbol z^-$ ) as the packet propagates outward, thereby initiating RDT. We initialise the fluctuations as linearly polarised Alfvén waves ( $\delta \boldsymbol{u}=\delta \boldsymbol{B}/\sqrt {4 \pi \rho }$ , corresponding to $\boldsymbol{z}^-=0$ ) with $\delta \boldsymbol u$ and $\delta \boldsymbol B$ lying in the $(\boldsymbol k\times {{\bar {\boldsymbol B}}})$ direction (Squire, Chandran & Meyrand Reference Squire, Chandran and Meyrand2020; Johnston et al. Reference Johnston, Squire, Mallet and Meyrand2022). The initial set of Alfvén waves are generated from a sum of randomly phased waves with amplitude determined by the Gaussian energy spectrum

(3.6) \begin{equation} E(k_\parallel ,k_\perp ) \propto \exp \!\Biggl [-\frac {(k_\parallel - k_{\parallel ,0})^2 + (k_\perp - k_{\perp ,0})^2}{k_w^2}\Biggr ]. \end{equation}

Here, $k_{\parallel ,0}=\kappa _\parallel ({2\pi }/{L_x}),\; k_{\perp ,0}=\kappa _\perp ({2\pi }/{L_\perp })\; \text{and}\; k_w=({12}/{L_\perp })$ , which set the scale of the centre of Gaussian peak and the width of the peak respectively; $k_{\parallel ,0}$ and $k_{\perp ,0}$ lie respectively parallel and perpendicular to $\bar {\boldsymbol B}$ , maintaining anisotropy even with an initial tilted spiral field. As the expanding-box simulation begins, the modes rapidly adjust within approximately 0.5 Alfvén times, losing the variation in $|\boldsymbol B|$ caused by their initial linear polarisation. This natural relaxation reduces the magnetic compressibility $C_B^2$ and allows the fluctuations to become spherically polarised.

To choose simulation parameters and quantify regimes, we use the time scales and their ratios, as discussed in § 2. The expansion time is $\tau _{\textrm {exp}}^{-1} = \dot a / a$ , the linear time scale is $\tau _{\textrm {A}}^{-1} = k_\parallel v_{\textrm {A}}$ and the outer-scale nonlinear time is $\tau _{\textrm {nl}}^{-1} = k_{\perp } z^+$ . From these rates we form the dimensionless control parameters: (i) $\chi _{\textrm {exp,0}} = (k_{\perp ,0} \; z^+_{\textrm {rms0}}) / (\dot a/a)$ ; (ii) $\chi _{\textrm {A}} \doteq (k_\perp z^+)/(k_\parallel v_{\textrm {A}})$ , which compares the strength of nonlinear interactions with linear Alfvénic propagation and thus quantifies how strongly nonlinearity competes with Alfvénic de-correlation of $z^-$ fluctuations (Goldreich & Sridhar Reference Goldreich and Sridhar1995); and (iii) the ratio of box-scale Alfvén time to the expansion time $\varDelta _0 \doteq \chi _{\textrm {exp0}}/\chi _{\textrm {A0}} = v_{\textrm {A,0}}/\dot {a} L_x$ which is constant. For the parameter set used in this paper we fix $\varDelta _0 = 2$ (via $\dot {a}=0.5$ , $L_\parallel =1$ and $v_{\textrm {A}}=1$ ) and vary $\chi _{\textrm {exp,0}}$ roughly in the range $\approx 11\!-\!30$ . We explore two initial fluctuation amplitudes defined by $A\equiv { z^+_{\textrm {rms0}}}/{v_{\textrm {A,0}}}=1.0 \;\text{and}\; 0.5,$ where $z^+_{\textrm {rms0}}$ is the initial root-mean-square fluctuation amplitude (see table 1 for exact values). The moderate-resolution simulations with reduced amplitude ( $A_0 = 0.5$ ) are designated as A05 to distinguish them from the standard moderate-resolution runs with $A_0 = 1.0$ . Changing $ A$ alters the initial nonlinear time (and therefore $\chi _{\textrm {exp,0}}$ and $\chi _{\textrm {A,0}}$ ), so we fix $\chi _{\textrm {exp,0}}$ to isolate pure amplitude effects. Physically, larger $A$ produces larger fluctuations with larger $\delta B_\parallel$ to maintain their spherical polarisation, whereas smaller $A$ (initial conditions) remain closer to the reduced MHD limit, which is derived assuming $z^+/v_{\textrm {A}}\ll 1$ .

3.3. Numerical issues at large $a$

We note that a numerical instability appears in our PS runs once the local spiral angle becomes very large ( ${\gtrsim} 70^\circ$ ), which occurs at sufficiently large $a$ . The instability manifests as small-scale, speckle-like noise in the transverse fields and an abrupt spike in inward Elsässer energy. The underlying cause is uncertain, but a plausible explanation is the increasing obliquity of the mean field along with the extreme anisotropic deformation of the EBM domain (the initially elongated box becoming a very flattened pancake), which under-resolves transverse modes and amplifies numerical error. This unfortunately stops us accessing any possible final balanced phase in PS runs. We therefore exclude affected intervals from our analysis. The issue is discussed in more detail in Appendix B.

5.1. Joint evolution of cross-helicity and residual energy

To quantify turbulence imbalance and the magnetic/kinetic partition, we track the normalised residual energy (4.1) and cross-helicity (4.2). We also define the alignment parameter $\sigma _\theta$

(5.1) \begin{equation} \sigma _\theta = \frac {\langle \boldsymbol{z}^+ \boldsymbol{\cdot }\boldsymbol{z}^- \rangle } {\sqrt {\langle |\boldsymbol{z}^+|^2 \rangle \, \langle |\boldsymbol{z}^-|^2 \rangle }} = \frac {\sigma _r}{\sqrt {1-\sigma _c^2}}\,, \end{equation}

such that fully aligned fluctuations with $\sigma _\theta =1$ lie on the circle $\sigma ^2_r + \sigma ^2_c = 1$ . Observations show fluctuations clustered near the lower-right quadrant of the unit circle, shifting from $\sigma _c\sim 1,\ \sigma _r\sim 0$ close to the Sun to $\sigma _c\sim 0,\ \sigma _r\sim -1$ at larger radii (Bruno et al. Reference Bruno, D’Amicis, Bavassano, Carbone and Sorriso-Valvo2007; Bruno & Carbone Reference Bruno and Carbone2013; Wicks et al. Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013). Figure 11 illustrates the trajectory of $(\sigma _c,\sigma _r)$ as the solar wind expands (parametrised by the expansion factor $a$ ), for both the pure radial field case and the PS case.

Figure 11.

Parametric evolution of cross-helicity $\sigma _{c}$ and residual energy $\sigma _{r}$ as a function of expansion factor $a$ . Panel (a) shows results from the HR- $\varPhi _0=0^\circ ,4^\circ$ simulations: circles correspond to the purely radial run and squares to the PS run. Panel (b) shows the A05- $\varPhi _0=0^\circ ,10^\circ$ simulations for both the radial and the higher initial Parker-angle case ( $\varPhi _0=10^\circ$ ). Coloured points represent $a = 1$ (dark blue), $a \sim 10$ (orange) and $a \sim 30$ (yellow). The black circle indicates the condition $\sigma _c^2 + \sigma _r^2 = 1$ . Inset of (a) shows the alignment parameter $\sigma _\theta$ versus $a$ for the HR– $\varPhi _0=0^\circ$ (blue) and HR– $\varPhi _0=4^\circ$ (orange) simulations. The radial run exhibits progressively stronger anti-alignment (more negative $\sigma _\theta$ ) with $a$ , while the PS case flattens, indicating saturation of the alignment. In (b) for the PS case, points with $a\gt 20$ are identified as affected by numerical instability: these points are plotted for completeness and the region they occupy is highlighted (highlighted, details in Appendix B).

In both cases, the system begins highly imbalanced ( $\sigma _c\approx 1$ , i.e. almost pure $z^+$ ) with $\sigma _r\approx 0$ . As $a$ increases, the curve stays close to the edge of the unit circle meaning $\sigma _r$ becomes increasingly negative. Physically, outward $\boldsymbol z^+$ fluctuations are reflected by the inhomogeneous wind into inward $\boldsymbol z^-$ that are nearly anti-aligned ( $\boldsymbol{z}^-\approx -\boldsymbol{z}^+$ ), so this naturally generates turbulence with $\langle \boldsymbol z^-\boldsymbol{\cdot }\boldsymbol z^+\rangle \lt 0$ (hence $\sigma _r\lt 0$ , (4.1)). This eventually produces magnetically dominated Alfvénic vortices (figure 3) as $\chi _{\textrm {exp}}$ decreases and $\sigma _c$ begins to decline toward zero (balanced turbulence), with $\sigma _r$ near its maximal negative values. The results are consistent with RMHD simulations of Meyrand et al. (Reference Meyrand, Squire, Mallet and Chandran2025) and in situ wind observations, which show transitions from imbalanced to nearly balanced turbulence with magnetic dominance beyond 1 AU (Bruno et al. Reference Bruno, D’Amicis, Bavassano, Carbone and Sorriso-Valvo2007).

The PS runs (squares) initially track the radial trend at small $a$ (the mean field is nearly radial there), but stay further from the circle edge at larger expansion. Even modest obliquity alters reflection and the polarisation of the reflected component: assuming an Alfvénic $\boldsymbol z^+$ predominantly perpendicular to $\bar {\boldsymbol B}$ , a reflected inward field $\boldsymbol z^-$ is no longer perfectly $-\boldsymbol z^+$ but instead a rotated combination of components due to the different sign of reflection in the $x$ direction (see (2.10)). This polarisation mismatch reduces the correlation $\langle \boldsymbol z^+\!\boldsymbol{\cdot }\boldsymbol z^-\rangle$ , so $\sigma _r$ remains less negative at a given $\sigma _c$ (which also decays more slowly with $a$ for the PS case, as shown in the inset of figure 5(a)). The inset of figure 11(a) shows $\sigma _\theta$ vs. $a$ for the HR- $\varPhi _{0}=0^\circ$ and HR- $\varPhi _{0}=4^\circ$ runs, quantifying the distance from the circle edge: HR- $\varPhi _{0}=0^\circ$ trends toward substantially stronger anti-alignment with increasing $a$ , whereas HR- $\varPhi _{0}=4^\circ$ saturates at a smaller $|\sigma _\theta |$ , which supports our view that the geometry-driven polarisation mismatch between the $yz$ an $x$ directions reduces alignment.

These results suggest an observational test: the joint evolution of $\sigma _c$ and $\sigma _r$ could be examined as a function of local PS angle. The spiral angle, defined by (2.1), varies with distance above the ecliptic (heliocentric latitude) at fixed radius. So, one prediction is that the high-latitude turbulence should hug the edge of circle more closely than turbulence in the ecliptic.

5.2. Switchbacks

Magnetic switchbacks (SBs) are sharp, transient reversals of the heliospheric magnetic-field direction that occur with little change in field magnitude. First detected by Helios and later observed in abundance by PSP, they are now recognised as a ubiquitous feature of the near-Sun solar wind (Bale et al. Reference Bale2019; Kasper et al. Reference Kasper2019; Raouafi et al. Reference Raouafi2023). Switchbacks are typically highly Alfvénic, their magnetic and velocity fields remain strongly correlated, and $|\boldsymbol B|$ is nearly constant throughout the reversal (Bale et al. Reference Bale2019).

Figure 12 illustrates representative fly-throughs at different stages of the evolution for both the radial and PS runs. In both geometries, $|\boldsymbol B|$ remains nearly constant while the field components undergo sharp directional swings in $B_\parallel$ , accompanied by correlated deflections in $B_T$ and $B_N$ . These signatures are consistent with large-amplitude, spherically polarised Alfvénic fluctuations as described in the past studies with similar set-ups (Squire et al. Reference Squire, Chandran and Meyrand2020; Shoda, Chandran & Cranmer Reference Shoda, Chandran and Cranmer2021; Johnston et al. Reference Johnston, Squire, Mallet and Meyrand2022). The emergence of nearly constant- $|\boldsymbol B|$ fluctuations reflects the nonlinear evolution’s tendency to favour spherically polarised, Alfvénic states that preserve field magnitude (Barnes & Hollweg Reference Barnes and Hollweg1974). This is likely because any perturbation that changes $|\boldsymbol B|$ excites compressive (fast or slow magnetosonic) modes with associated density and pressure variations; these compressive components steepen into shocks and dissipate more rapidly than the incompressible Alfvénic part that maintains constant field magnitude.

When $\delta B_\parallel$ reaches amplitudes comparable to the background, the local field direction can flip; the near constancy of $|\boldsymbol B|$ couples these radial reversals to transverse deflections, so each reversal appears as a coordinated multi-component swing. These rotations become noticeably sharper with expansion – compare the steeper swings between figures 12(a) and 12(c) – a natural consequence of the growing fluctuation amplitude. Furthermore, rotations in the PS runs are generally sharper than in the radial runs, as predicted by Squire et al. (Reference Squire, Johnston, Mallet and Meyrand2022). At late times, the radial case exhibits the loss of coordinated, constant- $|\boldsymbol{B}|$ rotations and increasingly irregular field-strength variations. These variations grow with distance, indicating that the magnetic field becomes less Alfvénic and more structurally complex as the cascade weakens. This late-time change broadly consistent with recent observational reports of enhanced field-strength variability with increasing heliocentric distance (Horbury et al. Reference Raouafi2023). In contrast, the PS configuration retains a higher degree of Alfvénicity at comparable expansion factors, maintaining near-constant $|\boldsymbol{B}|$ through large directional rotations and better preserving its nearly spherically polarised character.

Figure 12.

Fly throughs of the magnetic field for the high-resolution radial (HR– $\varPhi _0=0^\circ$ a, c, d) and PS (HR– $\varPhi _0=4^\circ$ , b, d, f) simulations at three expansion factors along the direction (1, 0.707, 0.392). Top row: (a) radial and (b) PS at $a=3.5$ ( $\varPhi \sim -13.75^\circ$ ); middle row: (c) radial and (d) PS at $a=11.5$ with $\varPhi \sim 38.8^\circ$ ; bottom row: (e) radial and (f) PS at $a=20$ at $\varPhi \sim -54.5^\circ$ . In each panel the total field strength $|\boldsymbol{B}|$ is shown in black; $B_\parallel$ , $B_T$ and $B_N$ are shown in blue, purple and orange, respectively. Each component is plotted normalised to $\bar {B}$ (vertical axis $B_i/\bar B$ ) versus normalised distance $l/(aL_\perp )$ , emphasising how component behaviour and $|\boldsymbol{B}|$ evolve with expansion.

We examine the SB fraction in our simulations, for which we define a SB as a region where the magnetic field $\boldsymbol B$ deviates from the local mean field $\bar {\boldsymbol B}$ by more than a specified threshold angle. This deviation is quantified by the normalised deflection parameter

(5.2) \begin{equation} z = \frac {1}{2}(1-\cos {\vartheta _z}), \end{equation}

where the deflection angle $\vartheta _z$ is given by

(5.3) \begin{equation} \cos {\vartheta _z} = \frac {\boldsymbol B \boldsymbol{\cdot }{\bar {\boldsymbol B}}}{\lvert \boldsymbol B\rvert \lvert {\bar {\boldsymbol B}}\rvert }. \end{equation}

Here, $z = 0$ corresponds to a field perfectly aligned with the background, and $z = 1$ to antiparallel configuration. We analyse regions corresponding to deflection thresholds $z_{\textrm {th}} = 0.125,\;0.25,\;0.375,\;0.5,\;0.625,\;0.75$ . We find that PS runs produce a larger fraction of strong directional reversals ( $z\gtrsim 0.5$ ), but a smaller fraction of weak deflections, than the corresponding radial runs (figure 13), in agreement with the three-dimensional simulations of Johnston et al. (Reference Johnston, Squire, Mallet and Meyrand2022) and analytic explanation of Squire et al. (Reference Squire, Johnston, Mallet and Meyrand2022). In the PS case, the SB fraction grows until a near spiral angle of $45^\circ$ and then declines. This decline coincides with the change in the evolution of $v_{\textrm {A}}$ between the two geometries discussed in § 2. The right panel of figure 13 shows the evolution of the normalised amplitude $z^+/v_{\textrm {A}}$ for all the simulation sets which are compared with the WKB expectations for a radial ( $z^+/v_{\textrm {A}} \propto a^{1/2}$ ) or azimuthal ( $z^+/v_{\textrm {A}} \propto a^{-1/2}$ ) background field. In the radial case, the amplitude rise slowly with increasing $a$ , whereas the PS runs show a different behaviour: they remain approximately flat or grow slowly at small $a$ and then decline at large $a$ .

Figure 13.

Evolution of SB fraction ( $f_{z} \geqslant z_{\textrm {th}}$ ) in MR- $\varPhi _0= 0^\circ$ (solid) and MR- $\varPhi _0=2^\circ$ (dashed lines). Both runs produce SBs, but the PS case exhibits systematically larger SB fractions across effectively all $a$ and a stronger growth of large-angle deflections with $a$ up to the point where fluctuation amplitudes decline; the downturn in ( $f_{z} \geqslant z_{\textrm {th}}$ ) at the largest $a$ is caused by the overall decrease in fluctuation amplitude. Right panel: the amplitude $z^+/v_{\textrm {A}}$ versus $a$ for all the simulations listed in the table 1. Dotted black lines show the WKB expectation for a radial field ( $z^+/v_{\textrm {A}} \propto a^{1/2}$ ) or azimuthally dominated field ( $z^+/v_{\textrm {A}} \propto a^{-1/2}$ ).

5.3. Compressibility

As noted in § 5.2, spherically polarised Alfvénic fluctuations are observed in the solar wind. To quantify the degree of compressibility in the evolving turbulence, we therefore measure the compressive fraction of the fluctuations. We use the diagnostics

(5.4) \begin{gather} {C_B = \sqrt {\frac {\langle (\delta |B|)^2\rangle }{\langle |\delta \boldsymbol{B}|^2\rangle }}} ,\quad \frac {\delta \rho _{\mathrm{rms}}}{\langle \rho \rangle } = \sqrt {\left \langle \left ( \frac {\rho - \langle \rho \rangle }{\langle \rho \rangle }\right )^2\right \rangle }. \end{gather}

The value of $C_B$ is nearly zero for pure transverse waves and grows when compressive fluctuations are present (Chen et al. Reference Chen2021; Shoda et al. Reference Shoda, Chandran and Cranmer2021) and $\delta \rho _{\textrm {rms}}/\langle \rho \rangle$ measures the fractional amplitude of density variations relative to the mean density. Taken together, $C_B$ and $\delta \rho _{\textrm {rms}}$ measure the prevalence of compressive fluctuations in the system. Figure 14 shows the evolution of $C_B$ and $\delta \rho _{\textrm {rms}}/\langle \rho \rangle$ as functions of the expansion factor $a$ . All runs start with higher $C_B$ because of the initialisation of simulations with large-amplitude linearly polarised waves. As expansion and nonlinear mode coupling proceed, compressions in $|\boldsymbol B|$ rapidly dissipated within $\leqslant 1 \tau _{\textrm {A}}$ , driving the system toward spherical polarisation and an effectively constant $|\boldsymbol B|$ ; correspondingly, $C_B$ drops significantly at early times.

Figure 14.

Evolution of (a) magnetic compressibility $C_B$ , and (b) density fluctuation amplitude $\delta \rho _{\textrm {rms}}/\langle \rho \rangle$ for the MR $-\varPhi _0=0^\circ ,2^\circ ,5^\circ$ simulations. All cases start with relatively large $C_B$ because of the linearly polarised initial conditions.

These diagnostics show no clear evidence that the PS geometry directly drives stronger compressive activity. Rather, PS runs retain higher normalised cross-helicity, (figure 5) producing sustained, near-constant- $|\boldsymbol B|$ rotations (Alfvénic SBs; figure 12). Compressive power increases as the turbulence becomes less Alfvénic – i.e. at lower $\sigma _c$ (or, equivalently, at larger $\sigma _r$ ) – as the relative contribution of $|\boldsymbol B|$ -varying compressive fluctuations grows (see figure 12(e)). Therefore, as the system evolves from a strongly imbalanced state toward a more balanced regime it loses its spherical polarisation and magnetic compressibility increases. Our simulations are consistent with this interpretation, with all cases similar at earlier times then differing later has $\sigma _c$ drops in the radial case, but remains higher for the PS ones (see figure 5).

We note an important caveat: the simulations presented here employ an locally isothermal equation of state, which does not capture the solar wind’s full thermodynamics. A more realistic equation of state should be included for direct quantitative comparison with spacecraft density signatures and compressive energetics.