1. Introduction
Alfvénic turbulence consists of non-compressive fluctuations in the velocity
$\boldsymbol{v}$
and magnetic field
$\boldsymbol{B}$
that have comparable energy densities and is believed to play a key role in generating the solar wind. The Sun launches Alfvén waves that travel outward and undergo partial reflection due to the radial variation in the Alfvén speed
$v_{\mathrm { A}}$
. The reflected waves then interact nonlinearly with the waves still traveling outward, causing fluctuation energy to cascade from large scales to small scales, where the fluctuations dissipate, heating and accelerating the plasma. Understanding the details of solar-wind turbulence has several important implications. First, it helps us predict how much mass the Sun loses through the solar wind and the speed at which this wind travels. The location where this turbulence heats the plasma is critical: heating inside the sonic critical point increases the mass loss rate, while heating farther away increases the wind’s final speed (Leer & Holzer Reference Leer and Holzer1980). Therefore, to determine the speed and mass flux of the solar wind, we need to understand how quickly this turbulence dissipates. Additionally, the inertial-range power spectrum of the turbulence, which determines the strength and anisotropy of small-scale fluctuations, is important to understand. These fluctuations are crucial because they control the amount of heating through mechanisms like stochastic heating and cyclotron heating (see, e.g. Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022). Finally, the degree of anisotropy in these fluctuations directly affects the efficiency of different dissipation mechanisms, which can control, for instance, the relative heating rate of different species (Johnston et al. Reference Johnston, Squire and Meyrand2024).
One of the most effective tools we have to study this turbulence is direct numerical simulations. A widely used approach involves flux-tube simulations, where we model a narrow magnetic flux tube centered on a radial magnetic-field line and track the turbulence as it evolves along the tube (Dmitruk & Matthaeus Reference Dmitruk and Matthaeus2003; van Ballegooijen et al. Reference van Ballegooijen, Asgari-Targhi, Cranmer and DeLuca2011; Perez & Chandran Reference Perez and Chandran2013; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2016, Reference van Ballegooijen and Asgari-Targhi2017; Chandran & Perez Reference Chandran and Perez2019). This method has shown great promise and has been extensively employed, but it presents a significant challenge: resolving the radial direction requires a vast number of grid points, leading to extremely high computational demands. A potential remedy is the so-called ‘expanding-box model’. The model simplifies the problem by focusing on a small, expanding parcel of plasma as it moves away from the Sun (Grappin et al. Reference Grappin, Velli and Mangeney1993; Dong et al. Reference Dong, Verdini and Grappin2014; Tenerani & Velli Reference Tenerani and Velli2017; Montagud-Camps et al. Reference Montagud-Camps, Grappin and Verdini2018; Shi et al. Reference Shi, Velli, Tenerani, Rappazzo and Réville2020; Squire et al. Reference Squire, Chandran and Meyrand2020; Johnston et al. Reference Johnston, Squire, Mallet and Meyrand2022; Grappin et al. Reference Grappin, Verdini and Müller2022; Meyrand et al. Reference Meyrand, Squire, Mallet and Chandran2025). The approach reduces computational demands while still capturing many essential characteristics of the turbulence. By concentrating on this localized region, the expanding-box model allows one to track how turbulence evolves as the plasma journeys through space, without resolving the full extent of the flux tube.
In this paper, we advance the expanding-box model by developing a generalized version tailored for reduced magnetohydrodynamics (RMHD). This new formulation allows the simulation box to either remain stationary in the Sun’s frame or move at any speed in the radial direction. Additionally, we have reformulated the equations using scalar potentials and Clebsch coordinates. This approach not only facilitates working with the flux-tube geometry but also enhances our understanding of the system while preserving the leading-order dynamics of the turbulent cascade. We believe that two types of simulations will be particularly valuable. First, the wave-frame box simulation, which follows the outward-propagating waves (
$z^+$
) that carry most of the energy, generalizes the standard expanding-box model to sub-Alfvénic regions close to the Sun and provides a complementary perspective on the development of turbulence. Second, the Eulerian box simulation at a fixed location in space will allow us to explore how turbulence evolves in a steady environment (Siggia & Patterson Reference Siggia and Patterson1978; Passot et al. Reference Passot, Sulem and Laveder2022).
The paper is structured as follows. In § 2, we present the mathematical foundation of our generalized expanding-box model. Section 3 details the scalar formulations of this model. Section 4 explores the linear dynamics, highlighting the significant effects of reflections compared with homogeneous RMHD. Finally, § 5 addresses the implications of our findings for understanding solar-wind turbulence and outlines potential directions for future research.
2. Governing equations
In the solar wind, the energetically dominant fluctuations are elongated along the magnetic field, to a degree that increases towards smaller scales (Horbury et al. Reference Horbury, Forman and Oughton2008; Sahraoui et al. Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010; Chen et al. Reference Chen, Mallet, Schekochihin, Horbury, Wicks and Bale2012). In the inertial range, in which the fluctuating magnetic field
$\delta \boldsymbol{B}$
is small compared with the root-mean-square magnetic-field strength
$B_{rms}$
, such anisotropic fluctuations are described by a version of RMHD (Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1974; Strauss Reference Strauss1976; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) that is modified to account for the outflow and inhomogeneity of the solar wind (e.g. Heinemann & Olbert Reference Heinemann and Olbert1980; Verdini & Velli Reference Verdini and Velli2007; Chandran & Hollweg Reference Chandran and Hollweg2009; van Ballegooijen et al. Reference van Ballegooijen, Asgari-Targhi, Cranmer and DeLuca2011; Perez & Chandran Reference Perez and Chandran2013; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2016). We refer to this set of equations as inhomogeneous RMHD, or IRMHD. In this section, we review the derivation of the IRMHD equations for solar-wind outflow within a narrow magnetic flux tube.
2.1. Narrow-flux-tube approximation
To properly describe an expanding plasma, we first need to design a magnetic flux tube. We assume it to be aligned along a radial magnetic-field line. Using spherical coordinates
$(r, \theta, \varphi )$
where
$r$
is the radial distance,
$\theta$
is the polar angle and
$\varphi$
is the azimuthal angle, we take
$\theta =0$
to coincide with the magnetic-field line upon which the flux tube centers and define
\begin{align} \bar {B}_0 \equiv B_0 (r, \theta = 0). \end{align}
We assume the background magnetic field
$B_0$
is independent of
$\varphi$
and has no
$\varphi$
component
\begin{align} B_{0 \varphi } = \frac {\partial B_{0 r}}{\partial \varphi } = \frac {\partial B_{0 \theta }}{\partial \varphi } = 0, \end{align}
and we require that
$\boldsymbol{B}_0$
is smooth around
$\theta =0$
(i.e.
$\nabla ^2 \boldsymbol{B}_0$
remains finite), which implies that
\begin{align} \left .\frac {\partial B_{0 r}}{\partial \theta }\right \vert _{\theta =0} = 0, \quad \text {and} \quad B_{0\theta }\Big \vert _{\theta =0} = 0. \end{align}
This results in
$B_{0 r} = \bar {B}_0 (r)+ O(\theta )^2$
. Given (2.2), the condition
$\boldsymbol {\nabla } \cdot \boldsymbol{B} = 0$
can be written as
\begin{align} \frac {1}{r^2} \frac {\partial }{\partial r} \left ( r^2 B_r\right ) + \frac {1}{r \sin \theta } \frac {\partial }{\partial \theta } \left ( B_\theta \sin \theta \right ) = 0. \end{align}
For a narrow flux tube,
$\theta \ll 1$
, and the polar component of the magnetic field can be expanded as
\begin{align} B_{0 \theta } = \sum _{n=1}^{\infty } c_n \theta ^n, \end{align}
where
$c_n$
are coefficients we need to determine, and the sum starts at
$n=1$
to satisfy (2.3). Upon substituting (2.5) into (2.4), we obtain
\begin{align} \frac {2 c_1}{r} + \frac {3 c_2 \theta }{r}= - \frac {1}{r^2} \frac {\partial }{\partial r} \left ( r^2 \bar {B}_0 \right ) + O(\theta )^2. \end{align}
From this, we find that
\begin{align} c_1 = - \frac {1}{r^2} \frac {\partial }{\partial r} \left ( r^2 \bar {B}_0 \right ), \quad \text {and} \quad c_2 = 0, \end{align}
indicating that the polar component of the magnetic field satisfies
\begin{align} B_{0 \theta } = - \frac {\theta }{2 r} \frac {\mathrm {d}}{\mathrm {d} r} \left ( r^2 \bar {B}_0 \right ) + O(\theta )^3. \end{align}
To summarize, we assume a background magnetic field
$\boldsymbol{B}_0 = \bar {B}_0 (r) \hat {\boldsymbol{e}}_r + B_{0\theta } \hat {\boldsymbol{e}}_\theta$
, with
$B_{0\theta }$
given by (2.8). This represents the narrow-flux-tube approximation, which has been extensively used to study waves and instabilities in flux tubes (see, e.g. Defouw Reference Defouw1976; Roberts & Webb Reference Roberts and Webb1978; Spruit Reference Spruit1981; van Ballegooijen et al. Reference van Ballegooijen, Asgari-Targhi, Cranmer and DeLuca2011). Following the work of Perez & Chandran (Reference Perez and Chandran2013), we also assume that the flux tube has a non-constant square cross-section and use Cartesian coordinates
$x$
and
$y$
, which are perpendicular to the central field line, such that
$\sqrt {x^2 + y^2} \ll r$
. This allows us to focus on the vicinity of the central field line. Finally, we neglect the Parker spiral effect, focusing our study on the inner heliosphere.
2.2. Inhomogeneous reduced magnetohydrodynamics
Our starting point is the ideal compressible magnetohydrodynamic equations
\begin{align} \frac {\partial \rho }{\partial t} &= - \boldsymbol {\nabla } \cdot \left ( \rho \boldsymbol{u} \right ), \end{align}
\begin{align} \rho \left ( \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \boldsymbol {\nabla } \boldsymbol{u}\right ) &= - \boldsymbol {\nabla } p_{\mathrm {tot}} + \frac {\boldsymbol{B} \cdot \boldsymbol {\nabla } \boldsymbol{B}}{4 \pi }, \end{align}
\begin{align} \frac {\partial \boldsymbol{B}}{\partial t} &= \boldsymbol {\nabla } \times \left ( \boldsymbol{u} \times \boldsymbol{B} \right ), \\[6pt] \nonumber \end{align}
where
$\rho$
is the mass density,
$\boldsymbol{u}$
is the velocity,
$\boldsymbol{B}$
is the magnetic field and
$p_{\mathrm {tot}} = P + B^2/8\pi$
is the the sum of the thermal and magnetic pressures, respectively. We set
\begin{align} \boldsymbol{u}(\boldsymbol{r},t) &= U(\boldsymbol{r}) \hat {\boldsymbol{b}}_0 + \delta \boldsymbol{u}(\boldsymbol{r},t), \\[-6pt] \nonumber \end{align}
\begin{align} \boldsymbol{B}(\boldsymbol{r},t) &= B_0(\boldsymbol{r}) \hat {\boldsymbol{b}}_0 + \delta \boldsymbol{B}(\boldsymbol{r},t),\\[-6pt] \nonumber \end{align}
\begin{align} \rho (\boldsymbol{r},t) &= \rho _0(\boldsymbol{r}), \\[6pt] \nonumber \end{align}
where the background flow velocity
$\boldsymbol{U}$
is aligned with the background magnetic field
$\boldsymbol{B}_0$
along the direction
$\hat {\boldsymbol{b}}_0 \equiv \boldsymbol{B}_0 / \left \vert \boldsymbol{B}_0 \right \vert$
, and we have neglected density fluctuations. Since we are interested in the Alfvén-wave dynamics, we assume transverse and non-compressive fluctuations, i.e.
\begin{align} \delta \boldsymbol{u} \cdot \boldsymbol{B}_0 &= \delta \boldsymbol{B} \cdot \boldsymbol{B}_0 = 0, \\[-6pt] \nonumber \end{align}
\begin{align} \boldsymbol {\nabla } \cdot \delta \boldsymbol{u} &= \boldsymbol {\nabla } \cdot \delta \boldsymbol{B} = 0, \\[6pt] \nonumber \end{align}
and we take
$\rho$
,
$U$
and
$B_0$
to be steady-state solutions of (2.9
a) through (2.9
c). The Alfvén velocity and Elsasser (Reference Elsasser1950) variables are respectively given by
\begin{align} \boldsymbol{v}_{\mathrm { A}} (\boldsymbol{r}) = \frac {\boldsymbol{B}_0 (\boldsymbol{r}) }{\sqrt {4 \pi \rho (\boldsymbol{r})}} \quad \text {and} \quad \boldsymbol{z}^\pm (\boldsymbol{r}, t)= \delta \boldsymbol{u} (\boldsymbol{r}, t) \mp \delta \boldsymbol{b} (\boldsymbol{r}, t), \end{align}
where
$\delta \boldsymbol{b} (\boldsymbol{r}, t) \equiv \delta \boldsymbol{B} (\boldsymbol{r}, t)/ \sqrt {4 \pi \rho (\boldsymbol{r})}$
. We take
$\boldsymbol{B}_0$
to be directed away from the Sun, so that
$\boldsymbol{z}^+$
(
$\boldsymbol{z}^-$
) corresponds to outward-propagating (inward-propagating) fluctuations. Given these assumptions, the fluctuations are well described by the equations of IRMHD (Velli et al. Reference Velli, Grappin and Mangeney1989; Zhou & Matthaeus Reference Zhou and Matthaeus1989; Cranmer & van Ballegooijen Reference Cranmer and van Ballegooijen2005; Verdini & Velli Reference Verdini and Velli2007; Perez & Chandran Reference Perez and Chandran2013)
\begin{align} \frac {\partial \boldsymbol{z}^\pm }{\partial t} + \left ( U \pm v_{\mathrm { A}} \right ) \left ( \hat {\boldsymbol{b}}_0 \cdot \boldsymbol {\nabla } \boldsymbol{z}^\pm \right ) + \left ( U \mp v_{\mathrm { A}} \right ) \left ( \frac {\boldsymbol{z}^\pm }{4 H_\rho } - \frac {\boldsymbol{z}^\mp }{2 H_{\mathrm { A}}} \right ) = - \boldsymbol{z}^\mp \cdot \boldsymbol {\nabla } \boldsymbol{z}^\pm - \boldsymbol {\nabla } p_\star, \end{align}
where
$p_\star = p_{\mathrm {tot}}/ \rho$
and the
$-\boldsymbol {\nabla } p_\star$
term on the right-hand side of (2.13) enforces
$ \boldsymbol {\nabla } \cdot \boldsymbol{z}^\pm = 0$
. The parameters
$H_\rho (r)$
, and
$H_{\mathrm { A}}(r)$
, are, respectively, the scale heights of the background mass density and the Alfvén speed. We also define the scale height of the background magnetic field
$H_{\mathrm { B}}(r)$
, giving
\begin{align} \frac {1}{H_{\mathrm { B}}} \equiv - \frac {\hat {\boldsymbol{b}}_0 \cdot \boldsymbol {\nabla } B_0}{B_0} \gt 0, \quad \frac {1}{H_\rho } \equiv - \frac {\hat {\boldsymbol{b}}_0 \cdot \boldsymbol {\nabla } \rho }{\rho } \gt 0, \quad \frac {1}{H_{\mathrm { A}}} \equiv \frac {\hat {\boldsymbol{b}}_0 \cdot \boldsymbol {\nabla } v_{\mathrm { A}}}{v_{\mathrm { A}}} = \frac {1}{2 H_\rho } - \frac {1}{H_{\mathrm { B}}}. \end{align}
Equation (2.13) is nicely expressed in terms of the wave-action variables introduced by Heinemann & Olbert (Reference Heinemann and Olbert1980)
\begin{align} \boldsymbol{f}^\pm (\boldsymbol{r}, t) \equiv \frac {1 \pm \eta ^{1/2}}{\eta ^{1/4}} \boldsymbol{z}^\pm, \quad \text {where} \quad \eta (\boldsymbol{r}) \equiv \frac {\rho (\boldsymbol{r})}{\rho _{\mathrm { A}}}, \end{align}
where
$\rho _{\mathrm { A}}$
is the density at the Alfvén surface, which is defined as the surface on which
$U = v_{\mathrm { A}}$
.Footnote
1
Upon substituting (2.14) and (2.15) into (2.13), we find that (Heinemann & Olbert Reference Heinemann and Olbert1980; Chandran & Hollweg Reference Chandran and Hollweg2009; Chandran & Perez Reference Chandran and Perez2019)
\begin{align} \frac {\partial \boldsymbol{f}^\pm }{\partial t} + \left ( U \pm v_{\mathrm { A}} \right ) \left ( \hat {\boldsymbol{b}}_0 \cdot \boldsymbol {\nabla } \boldsymbol{f}^\pm - \frac {\boldsymbol{f}^\mp }{2 H_{\mathrm { A}}} \right ) = - \boldsymbol{z}^\mp \cdot \boldsymbol {\nabla } \boldsymbol{f}^\pm - \frac {1 \pm \eta ^{1/2}}{\eta ^{1/4}} \boldsymbol {\nabla } p_\star . \end{align}
Thanks to the conservation of magnetic flux and mass,
$\rho U / B_0 =$
constant, and
$\eta = \mathcal {M}_{\mathrm { A}}^{-2}$
, where
$\mathcal {M}_{\mathrm { A}} \equiv U/v_{\mathrm { A}}$
is the Alfvénic Mach number. In the WKB (Wentzel—Kramers—Brillouin) limit, and in the absence of nonlinear interactions between counter-propagating waves, the average of
$\left (\boldsymbol{f}^+\right )^2$
over a wave period is independent of
$r$
(Heinemann & Olbert Reference Heinemann and Olbert1980), and hence, because of (2.15), the average of
$\left (\boldsymbol{z}^+\right )^2$
over a wave period peaks at the Alfvén critical point, where
$\mathcal {M}_{\mathrm { A}}=1$
.
2.3. Clebsch coordinates
To facilitate the design of future numerical simulations, we introduce Clebsch (Reference Clebsch1871) coordinates
$\tilde {x}$
and
$\tilde {y}$
that are constant along the background magnetic-field lines. In particular, we define
\begin{align} (\tilde {x}, \tilde {y}, \tilde {z})=(\alpha x, \alpha y, z), \end{align}
where
$\alpha ^2 \equiv \bar {B}_0/B_{\text {ref}}$
, and
$B_{\text {ref}}$
is a value of reference. The function
$\alpha (r)$
is the factor by which the separation between two field lines increases between radius
$r_{\mathrm { ref}}$
and some larger radius
$r$
. These coordinates satisfy the relation
\begin{align} \frac {\boldsymbol{B}_0}{B_{\text {ref}}} = \boldsymbol {\nabla } \tilde {x} \times \boldsymbol {\nabla } \tilde {y} + O(\theta ^2), \end{align}
from which it follows that
$\boldsymbol{B}_0 \cdot \boldsymbol {\nabla } \tilde {x} = \boldsymbol{B}_0 \cdot \boldsymbol {\nabla } \tilde {y} = 0$
. In the
$(\tilde {x}, \tilde {y}, \tilde {z})$
coordinate system, (2.13) and (2.16) transform into
\begin{align} \frac {\partial \boldsymbol{z}^\pm }{\partial t} + \left ( U \pm v_{\mathrm { A}} \right ) \frac {\partial \boldsymbol{z}^\pm }{\partial \tilde {z}} + \left ( U \mp v_{\mathrm { A}} \right ) \left ( \frac {\boldsymbol{z}^\pm }{4 H_\rho } - \frac {\boldsymbol{z}^\mp }{2 H_{\mathrm { A}}} \right ) = - \alpha \left ( \boldsymbol{z}^\mp \cdot \tilde {\boldsymbol {\nabla }} \boldsymbol{z}^\pm + \tilde {\boldsymbol {\nabla }}_\perp p_\star \right ), \end{align}
and
\begin{align} \frac {\partial \boldsymbol{f}^\pm }{\partial t} + \left ( U \pm v_{\mathrm { A}} \right ) \left ( \frac {\partial \boldsymbol{f}^\pm }{\partial \tilde {z}} - \frac {\boldsymbol{f}^\mp }{2 H_{\mathrm { A}}} \right ) = - \alpha \left ( \boldsymbol{z}^\mp \cdot \tilde {\boldsymbol {\nabla }} \boldsymbol{f}^\pm + \frac {1 \pm \eta ^{1/2}}{\eta ^{1/4}} \tilde {\boldsymbol {\nabla }}_\perp p_\star \right ), \end{align}
respectively, where
$\tilde {\boldsymbol {\nabla }}$
is the gradient operator in the
$(\tilde {x}, \tilde {y}, \tilde {z})$
coordinate system. In this coordinate system, the background magnetic-field lines are parallel to each other, as illustrated in figure 1. The parameter
$\alpha$
represents the separation distance between two field lines and adjusts the strength of the nonlinearities in the system. For the sake of readability, we will henceforth employ the Clebsch coordinates exclusively, and omit the tilde notation.
Figure 1. (a) A cartoon of a magnetic flux tube expanding from the Sun. The dark arrows represent magnetic-field lines, and the blue arrow shows the central field line of the flux tube. (b) A model of a narrow magnetic flux tube. (c) The same narrow magnetic flux tube after applying Clebsch coordinates. (d) Cartoon of the Eulerian box simulation: the box remains fixed while the plasma flows through it. (e) Cartoon of the moving box simulation: the box moves along with the plasma at a specified velocity, expanding with the field lines.
2.4. Transformation to a frame moving radially away from the Sun
In this section, we consider a reference frame that moves away from the Sun along the background magnetic-field lines at some arbitrary velocity
$U_{\mathrm { fr}}(z)$
. The position
$z(t)$
of a point that starts at
$z_0$
at time
$t_0$
and moves outward along
$\boldsymbol{B}_0$
at speed
$U_{\mathrm { fr}}(z)$
satisfies the equation
\begin{align} \frac {\partial }{\partial t} z \!\left (t \, | \, z_0, t_0 \right ) = U_{\mathrm {fr}} \left [ z \left (t \, | \, z_0, t_0 \right ) \right ]. \end{align}
In simpler words, this mapping relates the time elapsed since the initial moment
$t_0$
to the distance
$(z-z_0)$
traveled by the frame. With the aid of (2.21), we can transform (2.20) from the
$(x,y,z)$
coordinate system to the
$(x,y,z_0)$
coordinate system, obtaining
\begin{align} \frac {\partial \boldsymbol{f}^\pm }{\partial t} + \left ( U - U_{\mathrm {fr}} \pm v_{\mathrm { A}} \right ) \frac {U_{\mathrm {fr},0}}{U_{\mathrm {fr}}} \frac {\partial \boldsymbol{f}^\pm }{\partial z_0} - \frac {U \pm v_{\mathrm { A}}}{2 H_{\mathrm { A}}} \boldsymbol{f}^\mp = - \alpha \left ( \boldsymbol{z}^\mp \cdot \boldsymbol {\nabla } \boldsymbol{f}^\pm + \frac {1 \pm \eta ^{1/2}}{\eta ^{1/4}} \boldsymbol {\nabla } p_\star \right ), \end{align}
where
$U_{\mathrm {fr}}$
is shorthand for
$U_{\mathrm {fr}}(z(z_0,t))$
and
$U_{\mathrm {fr,0}} = U_{\mathrm {fr}}(z_0)$
The interested reader can find the details of the transformation calculations in Appendix A.
Equation (2.22) has two main differences compared with (2.16). First, there is a correction to the advection term:
$U \partial _z \to (U - U_{\mathrm {fr}}) \partial _z$
. This is intuitive because for a fixed observer (
$U_{\mathrm {fr}} = 0$
), we recover the usual advection term seen in the Eulerian frame. As the frame speed
$U_{\mathrm {fr}}$
gets closer to the plasma rest frame, the advection term decreases. Second, we must account for the changes in the
$z$
-derivative due to the motion of the frame. In general, as
$U_{\mathrm { fr}}$
depends on
$z$
, the coordinate system will be stretched or compressed. This deformation is captured by the transformation rule
$\partial /\partial z \rightarrow (U_{\mathrm {fr,0}}/U_{\mathrm {fr}}) \partial /\partial z_0$
.
Equation (2.22) represents the first main result of this paper. The main advantage of this formulation is its flexibility in choosing the appropriate frame of reference. For instance, in the plasma rest frame and the super-Alfvénic limit (
$U_{\mathrm {fr}} = U_{\mathrm {fr},0} = U \simeq \text {constant}$
and
$v_{\mathrm { A}} \ll U$
), we recover the system studied by Grappin et al. (Reference Grappin, Velli and Mangeney1993) and Meyrand et al. (Reference Meyrand, Squire, Mallet and Chandran2025). To revert to the Eulerian frame, we simply set the frame velocity to zero (
$U_{\mathrm {fr}} = U_{\mathrm {fr},0} = 0$
).
However, not all frames are equally relevant. While flux-tube simulations offer the most detailed insights into wave behavior, they are computationally expensive, and simplified models can be more efficient. Then, to understand how waves launched by the Sun evolve and to match the steady-state solution that arises in a full flux-tube simulation, it is natural to adopt the
$z^+$
frame, in which
$U_{\mathrm { fr}} = U + v_{\mathrm { A}}$
, and to allow the
$z^+$
fluctuations to decay in this frame without forcing
$\boldsymbol{z}^+$
. The resulting mean square value of
$z^+$
in the box,
$\langle z^+_{\mathrm { box}}(t)^2\rangle$
, can be used to estimate the average mean-square value of
$z^+$
in the solar wind,
$\langle z^+_{\mathrm { sw}}(r)^2\rangle$
, by setting
$\langle z^+_{\mathrm { sw}}(r)^2\rangle = \langle z^+_{\mathrm { box}}(t(r))^2\rangle$
, where
$t(r)$
is the time at which the moving box is at heliocentric distance
$r$
.
We note that, in a given simulation, if the box moves more slowly than
$U + v_{\mathrm { A}}$
, it will take longer to reach a given heliospheric distance, leading to excessive wave cascading and energy dissipation compared with that which could occur in the flux tube. To correct this, energy would have to be injected into the box. On the other hand, if the frame moves faster than
$U + v_{\mathrm { A}}$
, energy dissipation needs to be artificially increased. In the special case of a fixed Eulerian frame, the system must be forced so that
$z^+_{\mathrm {rms}}$
fluctuates about a steady mean, as it does at fixed
$r$
in the solar wind.
The Eulerian frame should provide a powerful method for investigating local turbulent processes, allowing us to focus on the system after it has reached a steady state. In this regime, time evolution becomes irrelevant, and computational resources can be allocated to enhancing resolution and extending the inertial range without the burden of tracking the dynamics of different spatial locations as required in moving box simulations. It is likely important in such simulations to choose the radial extent of the box to be longer than the distance that
$z^-$
fluctuations can propagate during their cascade time scales to prevent the same
$z^-$
and
$z^+$
fluctuations from interacting mutliple times. This approach should also enable the system to reach a steady state after just a few eddy turnover times, maximizing computational efficiency.
3. Scalar formulation
When solving the IRMHD equations on a computer, one can reduce the amount of computer time required by first rewriting the equations in terms of scalar potentials. To this end, we define the Elsasser potentials
$\zeta ^\pm$
via the equation
\begin{align} \boldsymbol{z}^\pm (\boldsymbol{x}, t) = \frac {\hat {\boldsymbol{e}}_{z} \times \boldsymbol {\nabla } \zeta ^\pm (\boldsymbol{x}, t)}{\alpha }. \end{align}
Equation (3.1) implies that the velocity and magnetic-field fluctuations can be expressed as
\begin{align} \delta \boldsymbol{u} (\boldsymbol{x}, t) = \frac {\hat {\boldsymbol{e}}_{z} \times \boldsymbol {\nabla } \Phi (\boldsymbol{x}, t)}{\alpha }, \quad \text {and} \quad \delta \boldsymbol{b} (\boldsymbol{x}, t) = \frac {\hat {\boldsymbol{e}}_{z} \times \boldsymbol {\nabla } \Psi (\boldsymbol{x}, t)}{\alpha }, \end{align}
where
$\Phi$
and
$\Psi$
are two stream functions that are determined up through an arbitrary additive constant and that satisfy the equation
$\zeta ^\pm = \left (\Phi \mp \Psi \right )$
. Similarly, we introduce the wave-action potentials
$\zeta _{\mathrm { f}}^\pm$
, which satisfy the equation
\begin{align} \boldsymbol{f}^\pm (\boldsymbol{x}, t) = \hat {\boldsymbol{e}}_{z} \times \boldsymbol {\nabla } \zeta _{\mathrm { f}}^\pm (\boldsymbol{x}, t), \quad \text {with} \quad \zeta _{\mathrm { f}}^\pm \equiv \frac {1 \pm \eta ^{1/2}}{\eta ^{1/4}} \zeta ^\pm . \end{align}
Finally, we define the field-aligned Elsasser and wave-action vorticities,
$\Omega ^\pm$
and
$\Omega _{\mathrm { f}}^\pm$
, via the equations
\begin{align} \Omega ^\pm (\boldsymbol{x}, t) &= \hat {\boldsymbol{e}}_{z} \cdot \left ( \boldsymbol {\nabla } \times \boldsymbol{z}^\pm \right ) = \nabla _\perp ^2 \zeta ^\pm, \end{align}
\begin{align} \Omega _{\mathrm { f}}^\pm (\boldsymbol{x}, t) &= \hat {\boldsymbol{e}}_{z} \cdot \left ( \boldsymbol {\nabla } \times \boldsymbol{f}^\pm \right ) = \nabla _\perp ^2 \zeta_{\mathrm{ f}}^\pm, \\[6pt] \nonumber \end{align}
where
$\nabla _\perp ^2 \equiv \partial _{xx} + \partial _{yy}$
. Throughout this paper, we henceforth assume
$\nabla _\perp \gg \partial _z \gg H_i^{-1}$
, with
$i = \{A, B, \rho \}$
. Given the previous assumptions, taking the curl of (2.22), and keeping only the terms of the lowest order, we obtain the second main result of this paper
\begin{align} \left [ \frac {\partial }{\partial t} + \left ( U - U_{\mathrm {fr}} \pm v_{\mathrm { A}} \right ) \frac {U_{\mathrm {fr},0}}{U_{\mathrm {fr}}} \frac {\partial }{\partial z_0} \right ] \Omega _{\mathrm { f}}^\pm - \frac {U \pm v_{\mathrm { A}}}{2 H_{\mathrm { A}}} \Omega _{\mathrm { f}}^\mp = \left \{ \partial _j \zeta _{\mathrm { f}}^\pm, \partial _j \zeta ^\mp \right \} + \left \{ \Omega _{\mathrm { f}}^\pm, \zeta ^\mp \right \}, \end{align}
where we have employed the Einstein summation convention with the index
$j=\{x,y\}$
, and where the Poisson bracket
$\left \{f, g\right \} \equiv \hat {\boldsymbol{e}}_{z} \cdot \left (\boldsymbol {\nabla } f \times \boldsymbol {\nabla } g \right )$
for arbitrary functions
$f$
and
$g$
. (3.5) can be divided in three parts:
-
(i) The squared brackets on the left-hand side correspond to a linear advection of wave-action vorticities
$\Omega _{\mathrm { f}}^\pm$
at speed
$\left ( U - U_{\mathrm {fr}} \pm v_{\mathrm { A}} \right ) U_{\mathrm {fr},0} / U_{\mathrm {fr}}$
. -
(ii) The last term on the left-hand side describes non-WKB reflections, which act as a source for inward-propagating waves that are correlated to the original wave. This correlation leads to the phenomenon often referred to as anomalous coherence (see, e.g. Velli et al. Reference Velli, Grappin and Mangeney1989; Verdini et al. Reference Verdini, Velli and Buchlin2009; Perez & Chandran Reference Perez and Chandran2013; Chandran & Perez Reference Chandran and Perez2019). This term also disrupts the conservation of energy and cross helicity for the turbulent fluctuations. However, these conservation laws are restored when the energy and cross helicity of the background are included (Chandran et al. Reference Chandran, Perez, Verscharen, Klein and Mallet2015).
-
(iii) The nonlinear terms capture two key effects. First (by order of appearance), the vortex stretching effect, arising from interactions between opposite populations, influences the growth or decay of vorticities (Zhdankin et al. Reference Zhdankin, Boldyrev and Uzdensky2016). As highlighted by Schekochihin (Reference Schekochihin2022), the nonlinear term forces
$\Omega _{\mathrm { f}}^+$
and
$\Omega _{\mathrm { f}}^-$
equally but in opposite directions at every point in space and time, resulting in a negative correlation between the two vorticities. This interaction systematically promotes current sheets over shear layers. Second, nonlinear advection describes how each vorticity is transported by the field of the other.
We note that (3.5) becomes problematic near the Alfvén critical point
$r_{\mathrm { A}}$
, where
$\Omega _{\mathrm { f}}^-$
and
$\zeta _{\mathrm { f}}^-$
vanish, and where the value of
$\zeta ^-$
in the equation for
$\partial \Omega _{\mathrm { f}}^+/\partial t$
would need to be computed by dividing
$\zeta _{\mathrm { f}}^-$
by
$1-\eta ^{1/2}$
, which also vanishes at
$r=r_{\mathrm { A}}$
. To avoid the need to evaluate
$\zeta _{\mathrm { f}}^-/(1-\eta ^{1/2})$
in moving-box simulations that transit past
$r=r_{\mathrm { A}}$
, one could replace the equation for
$\partial \Omega _{\mathrm { f}}^-/\partial t$
by an equivalent equation for
$\partial \Omega ^-/\partial t$
.Footnote
2
4. Linear system
In this subsection, we consider a rectangular box at some initial position in the Sun’s frame. We take the radial extent
$\Delta r$
of the box to be much smaller than the distance over which the equilibrium properties change appreciably. Anticipating future numerical simulations, we make the approximation that the fluctuations in the box satisfy triply periodic boundary conditions. We further take
$H_{\mathrm { A}}$
to be a non-zero constant, but we treat
$v_{\mathrm { A}}$
and
$U$
as uniform within the box. Although a spatially uniform
$v_{\mathrm { A}}$
is inconsistent with a non-zero
$H_{\mathrm { A}}$
, the variations in
$v_{\mathrm { A}}$
that we are neglecting are
$\sim \Delta r / H_{\mathrm { A}} \ll 1$
, and, we conjecture, unimportant—similarly to Boussinesq’s approximation. In the general case, simple linear sinusoidal solutions to (3.5) cannot be obtained because
$U(z(t))$
, and
$v_{\mathrm { A}}(z(t))$
depends on time in a non-trivial way.Footnote
3
We therefore choose to focus this section on the Eulerian frame, in which the box remains at a fixed position in the Sun’s frame. We now consider a small perturbation of the wave-action potentials
$\zeta _{\mathrm { f}}^\pm = \zeta ^\pm _{\mathrm { f}0} \mathrm {e}^{ \mathrm {i} \left ( \boldsymbol{k} \cdot \boldsymbol{x} - \omega t \right )}$
. Upon substituting this perturbation into (3.5), we obtain the matrix equation
\begin{align} \left(\begin{array}{c@{\quad}c} \mathrm {i} \left ( \omega ^{\prime} - \omega _{\mathrm { A}} \right ) & (U + v_{\mathrm { A}})(2H_{\mathrm { A}})^{-1} \\[5pt] (U - v_{\mathrm { A}})(2H_{\mathrm { A}})^{-1} & \mathrm {i} \left ( \omega ^{\prime} + \omega _{\mathrm { A}} \right ) \end{array} \right) \begin{pmatrix} \zeta ^+_{\mathrm { f}0} \\[3pt] \zeta ^-_{\mathrm { f}0} \\ \end{pmatrix} = 0, \end{align}
with
$\omega ^{\prime} \equiv \omega - k_z U$
, and
$\omega_{\mathrm{A}} \equiv k_z v_{\mathrm { A}}$
. This linear set of equations supports waves forward- and backward-propagating modes of frequency
\begin{align} \omega _\pm ^{\prime}(k_z) = \pm v_{\mathrm { A}} \sqrt { k_z^2 + \frac {1 - \mathcal {M}_{\mathrm { A}}^2}{4 H_{\mathrm { A}}^2 }}. \end{align}
This equation generalizes the Alfvén-wave dispersion relation to allow for finite Alfvén-speed scale heights, which makes the waves dispersive. For non-zero values of
$k_\|$
, the second term under the square-root sign in (4.2) is a small correction to the first term as
$\Delta r \ll H_{\mathrm { A}}$
. However, when
$k_\|=0$
, that second term is the only term that survives, leading to purely growing or damped solutions when
$M_A \gt 1$
, and to low-frequency oscillatory solutions when
$M_A \lt 1$
.
We note that the eigenfunctions that diagonalize the linear system are
\begin{align} \Theta _k^\pm = \mp \frac {U-v_{\mathrm { A}}}{2 H_{\mathrm { A}}} \xi _k^+ + \mathrm {i} \left ( \omega _+^{\prime} \mp \omega _{\mathrm { A}} \right ) \xi _k^-. \end{align}
Unlike RMHD, in which the eigenmodes are
$\zeta ^\pm _k$
, here,
$\Theta _k^\pm$
combine both
$\xi ^+_{k}$
and
$\xi ^-_{k}$
. This distinction introduces an important new feature: the governing equations for the eigenfunctions
$\Theta _k^\pm$
inherently describe nonlinear interactions between both co-propagating and counter-propagating waves.
5. Discussion and conclusion
In this paper, we have developed a family of local approximations to the nonlinear equations derived by Heinemann & Olbert (Reference Heinemann and Olbert1980), including a box that moves radially outward at an arbitrary velocity. The use of Clebsch coordinates ensures independence from radial magnetic-field line spreading and greatly simplifies the system’s geometry, allowing us to derive (3.5), which is a relevant extension for both sub- and super-Alfvénic plasmas, regardless of whether the turbulence is balanced or imbalanced, and whether the medium is homogeneous or not.
There are several numerical approaches to solving (3.5), each with its advantages and limitations. Flux-tube simulations, for example, are the most comprehensive as they cover the full relevant range in
$r$
, making them ideal for capturing the full complexity of the turbulence (see e.g. Perez & Chandran Reference Perez and Chandran2013; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2016, Reference van Ballegooijen and Asgari-Targhi2017; Chandran & Perez Reference Chandran and Perez2019). However, they are incredibly expensive computationally. This is because they require an extremely large number of radial grid points, and also because Alfvén waves must travel from the Sun out to the Alfvén critical point and back a few times before the turbulence settles down into a statistical steady state (Perez & Chandran Reference Perez and Chandran2013). The Eulerian box, a small stationary box in the expanding solar wind, will allow a wide inertial range and is particularly useful for studying local turbulence dynamics and phenomena like cascade rates and anomalous coherence. The trade-off, however, is that it requires forcing to sustain the root mean square value of
$\boldsymbol{z}^+$
, denoted
$z^+_{rms}$
, at a chosen level, and
$z^+_{\mathrm { rms}}$
and the parallel box length must be large enough that
$\boldsymbol{z}^-$
fluctuations to a good approximation cascade and dissipate before crossing the box to avoid spurious interactions with their parent outward-propagating waves. Additionally, it lacks radial tracking, so while one can get insight into inertial-range dynamics, it misses out on how the turbulence evolves with heliospheric distance and thus cannot predict coronal/solar-wind heating rates. Finally, the moving box method generalizes the well-known expanding-box model to a non-constant solar wind’s speed, and strikes a balance between these two approaches. This approach corresponds to the accelerating expanding-box model derived by Tenerani & Velli (Reference Tenerani and Velli2017) in the specific limit of incompressible and strongly anisotropic fluctuations. However, our model introduces a key advantage: it offers the flexibility to select a reference frame that is best suited to the problem at hand. In contrast, the accelerating expanding-box model is inherently tied to the plasma bulk flow, which becomes less useful in the sub-Alfvénic wind, as discussed in § 2.4. Like the Eulerian box, the moving box method gives a large inertial range and can explore anomalous coherence, but it also captures the decay of turbulence as it propagates away from the Sun, in a way that is independent of the forcing mechanism—if the
$z^+$
frame is used. One of its main strengths is the ability to track turbulence radially as the simulation progresses. However, it comes with the drawback of limited statistics at any given radius, since the box moves along with the plasma flow. Each method offers its own unique insights, but the choice ultimately depends on the specific aspect of turbulence one is trying to study.
Finally, our focus has been on Alfvén waves. To extend this work we could cover the Alfvénic cascade down to electron scales, deriving an extension of finite Larmor radius MHD (Schekochihin et al. Reference Schekochihin, Kawazura and Barnes2019; Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021) to finite scale heights.
Acknowledgements
Editor Thierry Passot thanks the referees for their advice in evaluating this article.
Funding
This work was supported in part by NASA grant 80NSSC24K0171 and also NASA grant NNN06AA01C to the Parker Solar Probe FIELDS Experiment. J.S. acknowledges the support of the Royal Society Te Apārangi through Marsden-Fund grant MFP-UOO2221, as well as through the Rutherford Discovery Fellowship RDF-U001804.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Moving frame transformation
To understand the transformation to a frame moving at speed
$U_{\mathrm {fr}}(z)$
, we introduce new time and radial variables
$t^{\prime}$
and
$z_0$
, with
$t^\prime = t$
. The coordinates
$x$
and
$y$
remain unchanged and, unlike the time variable, are not coupled to the transformation, so we ignore them here for brevity. We start by defining the radial position
$z$
of a point initially located at position
$z_0$
at time
$t_0^{\prime}$
as a function of
$t^{\prime}$
, written as
$z \left (t^{\prime} \, | \, z_0, t_0^{\prime} \right )$
, assuming the point moves outward at speed
$U_{\mathrm {fr}}(z)$
. The evolution of
$z$
over time is given by
\begin{align} \frac {\partial }{\partial t^{\prime}} z \left (t^{\prime} \, | \, z_0, t_0^{\prime} \right ) = U_{\mathrm {fr}} \left [ z \left (t^{\prime} \, | \, z_0, t_0^{\prime} \right ) \right ]. \end{align}
This equation can be integrated to yield
\begin{align} z \left (t^{\prime} \, | \, z_0, t_0^{\prime} \right ) - z_0 = \int _0^{t^{\prime}} U_{\mathrm {fr}} \left [ z \left (t_1 \, | \, z_0, t_0^\prime \right ) \right ] \mathrm {d}t_1. \end{align}
Likewise, if the initial position of the point at time
$t_0^\prime$
is
$z_0 + \Delta$
, where
$\Delta$
is some small radial displacement, then the point’s position at time
$t^\prime$
satisfies the equation
\begin{align} z \left (t^{\prime} \, | \, z_0 + \Delta, t_0^{\prime} \right ) - z_0 - \Delta = \int _0^{t^{\prime}} U_{\mathrm {fr}} \left [ z \left (t_1 \, | \, z_0 + \Delta, t_0^\prime \right ) \right ] \mathrm {d}t_1. \end{align}
We define
$\delta t$
such that
$z \left (t_0^{\prime} + \delta t \, | \, z_0, t_0^{\prime} \right ) = z_0 + \Delta .$
That is,
$\delta t$
is the amount of time required for a point starting at
$z_0$
to reach position
$z_0 + \Delta$
. It then follows that
\begin{align} z \left (t^{\prime} + \delta t \, | \, z_0, t_0^{\prime} \right ) = z \left (t^{\prime} + \delta t \, | \, z_0 + \Delta, t_0^{\prime} + \delta t\right ) = z \left (t^{\prime} \, | \, z_0 + \Delta, t_0^{\prime} \right ). \end{align}
The first equality in (A.4) results from the fact that the point that starts at
$(z_0, t^\prime _0)$
moves through space–time coordinates
$(z_0 + \Delta, t^\prime _0 + \delta t)$
on its way to final position
$z \left (t^{\prime} + \delta t \, | \, z_0, t_0^{\prime} \right )$
. The second equality in (A.4) follows from the fact that
$U_{\mathrm { fr}}(z)$
depends on
$z$
but not
$t^\prime$
. When
$\delta t$
and
$\Delta$
are infinitesimal
\begin{align} \delta t = \frac {\Delta }{U_{\mathrm {fr}}(z_0)}. \end{align}
Using (A.4) to rewrite the right-hand side of (A.3) and then subtracting (A.2) from (A.3), we obtain
\begin{align} z(t^\prime \, |\, z_0 + \Delta, t_0^{\prime}) - z(t^\prime \,|\,z_0, t_0^{\prime}) - \Delta & = \int _{t^\prime }^{t^\prime + \delta t} U_{\mathrm {fr}} \left [ z \left (t_1 \, | \, z_0, t_0^\prime \right ) \right ] \mathrm {d}t_1 \nonumber\\ & \quad - \int _0^{\delta t} U_{\mathrm {fr}} \left [ z \left (t_1 \, | \, z_0, t_0^\prime \right ) \right ] \mathrm {d}t_1. \end{align}
Noting that the final terms on the left- and right-hand sides of (A.6) cancel and dividing (A.6) by
$\Delta$
, we arrive at
\begin{align} \frac {z(t^\prime \, | \, z_0 + \Delta, t_0^{\prime}) - z(t^\prime \, | \, z_0, t_0^{\prime})}{\Delta } = \frac {\delta t}{\Delta } U_{\mathrm {fr}} \left [ z \left (t^\prime \, | \, z_0, t_0^\prime \right ) \right ]. \end{align}
Finally, upon taking the limit
$\Delta \rightarrow 0$
, we find that
\begin{align} \left .\frac {\partial z}{\partial z_0}\right |_{t^{\prime}} = \frac {U_{\mathrm {fr}}(z)}{U_{\mathrm {fr}}(z_0)} = \frac {U_{\mathrm {fr}}}{U_{\mathrm {fr},0}}. \end{align}
Equation (A.8) describes how the coordinate system is stretched or compressed by the radial variations in
$U_{\mathrm { fr}}(z)$
. The Jacobian matrix associated with this change of variables is
\begin{align} J \equiv \frac {\partial \left ( z, t \right ) }{\partial \left ( z_0, t^{\prime} \right )} = \left(\begin{array}{l@{\quad}l} \left .\frac {\partial z}{\partial z_0}\right |_{t^{\prime}} & \left .\frac {\partial z}{\partial t^{\prime}}\right |_{z_0}\\[6pt] \left .\frac {\partial t}{\partial z_0}\right |_{t^{\prime}} & \left .\frac {\partial t}{\partial t^{\prime}}\right |_{z_0} \end{array} \right) = \left(\begin{array}{c@{\quad}c} \frac {U_{\mathrm {fr}}}{U_{\mathrm {fr},0}} & U_{\mathrm {fr}}\\ 0 & 1 \end{array}\right) .\end{align}
It follows from the chain rule that the Jacobian matrix of the inverse transformation is the inverse of the matrix in (A.9). Hence,
\begin{align} J^{-1} = \frac {\partial \left ( z_0, t^{\prime} \right )}{\partial \left ( z, t \right ) } = \left(\begin{array}{c@{\quad}c} \left .\frac {\partial z_0}{\partial z}\right |_{t} & \left .\frac {\partial z_0}{\partial t}\right |_{z}\\[6pt] \left .\frac {\partial t^{\prime}}{\partial z}\right |_{t} & \left .\frac {\partial t^{\prime}}{\partial t}\right |_{z} \end{array} \right) = \left(\begin{array}{l@{\quad}l} \frac {U_{\mathrm {fr},0}}{U_{\mathrm {fr}}} & -U_{\mathrm {fr}}\\ 0 & 1 \end{array} \right), \end{align}
where the final equality follows from taking the inverse of the
$2\times 2$
matrix on the right-hand side of (A.9). Equation (A.10) provides the partial derivatives that we will need in order to transform (2.20) to the
$(z_0, t^\prime )$
coordinate system. In particular, (A.10) and the chain rule enable us to write
\begin{align} \frac {\partial \boldsymbol{f}^\pm }{\partial t} &= \left .\frac {\partial \boldsymbol{f}^\pm }{\partial t^{\prime}}\right |_{z_0} - U_{\mathrm {fr},0} \frac {\partial \boldsymbol{f}^\pm }{\partial z_0}, \\[-6pt] \nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{f}^\pm }{\partial z} &= \left .\frac {\partial t^{\prime}}{\partial z}\right |_{t} \left .\frac {\partial \boldsymbol{f}^\pm }{\partial t^{\prime}}\right |_{z_0} + \left .\frac {\partial z_0}{\partial z}\right |_{t} \left .\frac {\partial \boldsymbol{f}^\pm }{\partial z_0}\right |_{t^{\prime}} = \frac {U_{\mathrm {fr},0}}{U_{\mathrm {fr}}} \frac {\partial \boldsymbol{f}^\pm }{\partial z_0}. \\[6pt] \nonumber \end{align}
Upon substituting these expressions into (2.20) and omitting the prime symbol on
$t^\prime$
to make the notation more compact, we obtain (2.22).
Open access
