# Higher-Order Analysis of Three-Dimensional Anisotropy in Imbalanced Alfvénic
Turbulence

Nikos Sioulas Department of Earth, Planetary, and Space Sciences, University
of California, Los Angeles, CA, USA  Themistocles Zikopoulos Department of
Physics, Aristotle University of Thessaloniki, GR-52124 Thessaloniki, Greece
Chen Shi Department of Earth, Planetary, and Space Sciences, University of
California, Los Angeles, CA, USA  Marco Velli Department of Earth, Planetary,
and Space Sciences, University of California, Los Angeles, CA, USA  Trevor
Bowen Space Sciences Laboratory, University of California, Berkeley, CA
94720-7450, USA  Alfred Mallet Space Sciences Laboratory, University of
California, Berkeley, CA 94720-7450, USA  Luca Sorriso-Valvo Institute for
Plasma Science and Technology (ISTP), CNR, Bari, Italy  Space and Plasma
Physics, School of Electrical Engineering and Computer Science, KTH Royal
Institute of Technology, Stockholm, Sweden  Andrea Verdini Universitá di
Firenze, Dipartimento di Fisica e Astronomia, Firenze, Italy  B. D. G.
Chandran Space Science Center and Department of Physics, University of New
Hampshire, Durham, NH 03824, USA  Mihailo M. Martinović Lunar and Planetary
Laboratory, University of Arizona, Tucson, AZ 85721, USA  S. S. Cerri
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, Bd de l’Observatoire, CS 34229, 06304 Nice5cedex 4, France  Nooshin
Davis Space Science Center and Department of Physics, University of New
Hampshire, Durham, NH 03824, USA  Corina Dunn Space Sciences Laboratory,
University of California, Berkeley, CA 94720-7450, USA

###### Abstract

In-situ observations of imbalanced solar wind turbulence are employed to
evaluate the core assumptions and higher-order scaling predictions of MHD
turbulence models grounded in the principles of “Critical Balance” (CB) and
“Scale-Dependent Dynamic Alignment” (SDDA).

Our results indicate that, at the energy injection scales, both outgoing and
ingoing Alfvénic fluctuations undergo a weak energy cascade,
$\chi_{\lambda}^{\pm}<1$, where
$\chi_{\lambda}^{\pm}\equiv\tau_{A}^{\pm}/\tau_{nl}^{\pm}$, the ratio of
linear to non-linear timescales. Simultaneously, a tightening of SDDA is
observed across this range. While the outgoing waves remain in a weak cascade
state throughout the inertial range, $\chi_{\lambda}^{+}\approx 0.2$, the
ingoing waves transition to a strong cascade, $\chi_{\lambda}^{-}>1$, at
$\lambda\approx 3\times 10^{4}d_{i}$. This transition, however, is accompanied
by spectral scalings that diverge from the expected canonical behavior marking
the shift from weak to strong turbulence - a discrepancy we attribute to the
effects of “anomalous coherence”. The domain canonically identified as the
inertial range consists of two distinct sub-inertial segments. At
$\lambda\gtrsim 100d_{i}$, the “average” eddy assumes a field-aligned tube
topology, with SDDA signatures being weak and largely restricted to the
highest amplitude fluctuations. The scaling exponents $\zeta_{n}$ of the
$n$-th order conditional structure functions, perpendicular to both the local
mean field and the fluctuation direction, conform to the analytical models of
Chandran et al. (2015) and Mallet & Schekochihin (2017), indicating
“multifractal” statistics and strong intermittency; the scaling in the
parallel and displacement (i.e., fluctuation direction) components is more
concave than theoretically predicted. We argue that the statistics of this
range could be contaminated by expansion effects. Below $\lambda\approx
100d_{i}$, eddies display increasing anisotropy, evolving into structures
resembling thin current sheets. Concurrently, $\zeta_{n}$ scales linearly with
order, signaling a transition towards “monofractal” statistics. At
$\lambda\approx 8d_{i}$, the increase in aspect ratio ceases, and the eddies
transition to a quasi-isotropic state. This shift might be a signature of the
tearing instability, potentially leading to reconnection of the thin current
sheets, or it could result from turbulent energy being channeled into an ion-
cyclotron wave spectrum, consistent with the “helicity barrier” effect.

Our analysis employs 5-point structure functions, which are shown to be more
effective than the traditional 2-point approach in accurately capturing the
steep scaling behaviors observed at smaller spatial scales.

Magnetohydrodynamics(MHD); Solar Wind; Plasmas; Turbulence; Waves

††software: Python (Van Rossum & Drake Jr, 1995), SciPy (Virtanen et al.,
2020), Pandas (McKinney et al., 2010), Matplotlib (Hunter, 2007), Pyspedas
(Angelopoulos et al., 2019), (Sioulas, 2023)

##  1 Introduction

Propelled by the internal dynamics of the Sun, the solar wind expands
spherically within the heliosphere, transporting a wide spectrum of magnetic
field and plasma fluctuations (Bruno & Carbone, 2013). Given that fluctuations
are predominantly weakly compressive (Tu & Marsch, 1995), their dynamics are
often analyzed within the framework of incompressible magnetohydrodynamic
(MHD) theory (Biskamp, 2003; Matthaeus & Velli, 2011). The Elsasser (1950)
form of the ideal incompressible MHD equations is given by

$\frac{\partial\delta\bm{z}^{\pm}}{\partial
t}\mp(\bm{V_{A}}\cdot\nabla)\delta\bm{z}^{\pm}=-(\delta\bm{z}^{\mp}\cdot\nabla)\delta\bm{z}^{\pm}-\frac{\nabla
p}{\rho_{0}},$ (1)

where the Elsässer variables,
$\delta\bm{z}^{\pm}=\delta\bm{v}\pm\delta\bm{b}/\sqrt{4\pi\rho}$, represent
eigenfunctions of Alfvén (1942) waves propagating (anti)parallel to the
background magnetic field ($\bm{B}_{0}$) at the Alfvén speed,
$\bm{V}_{a}=\bm{B}_{0}/\sqrt{4\pi\rho}$. The total energy $E_{t}=E^{+}+E^{-}$,
and cross-helicity $H_{c}=E^{+}-E^{-}$, expressed in terms of the energy
associated with fluctuations in $\bm{z}^{\pm}$,
$E^{\pm}=\langle|\delta\bm{z}^{\pm}|^{2}\rangle/4$ are ideal (i.e., with zero
viscosity and resistivity) invariants of the incompressible MHD equations.
Their ratio defines the normalized cross helicity $\sigma_{c}=H_{c}/E_{t}$.
When the energy fluxes, denoted as $\epsilon^{\pm}$, in wave packets traveling
in opposite directions differ, indicated by $\epsilon^{+}/\epsilon^{-}\neq 1$,
turbulence is termed as imbalanced, $\sigma_{c}\neq 0$. The nonlinearities, as
highlighted by the first term on the right-hand side of Equation 1, arise from
collisions between oscillation modes with opposite signs of cross-helicity
(i.e. between $\delta\bm{z}^{+}$ and $\delta\bm{z}^{-}$). These inertial, or
energy-conserving, interactions lead to the fragmentation of the wavepackets
and drive energies of $E^{+}$ and $E^{-}$ towards smaller perpendicular
scales, at which point energy is converted into heat through dissipation
(Iroshnikov, 1963; Kraichnan, 1965), henceforth referred to as IK. The
turbulent cascade is strongly anisotropic and results in the formation of
structures characterized by $\ell_{||}\gg\lambda$, with $\ell_{||}\sim
1/k_{||}$ and $\lambda\sim 1/k_{\perp}$ representing the correlation lengths
along and perpendicular to the background magnetic field111In assessing the
dominant scale of the background magnetic field, two methodologies are
predominantly employed, utilizing the global and local frames (Cho & Vishniac,
2000; Maron & Goldreich, 2001). For a detailed exploration of the consequences
associated with defining the magnetic field either globally, thereby
$\bm{B}_{0}$, or locally, denoted as $\bm{B}_{\ell}$, interested readers are
encouraged to refer to (Chen et al., 2011; Matthaeus et al., 2012; Gerick et
al., 2017)., respectively (Montgomery & Turner, 1981; Shebalin et al., 1983;
Grappin, 1986). The effectiveness of the collisions hinges on two critical
timescales: the time it takes for a turbulent perturbation of size $\lambda$
to break up nonlinearly $\tau_{nl}^{\pm}\sim\lambda/\delta z^{\mp}$ and the
linear propagation (or collision) time $\tau_{A}^{\pm}=\ell_{||}^{\pm}/V_{a}$.
The ratio of these timescales defines the nonlinearity parameter,
$\chi^{\pm}\equiv\tau_{A}^{\pm}/\tau_{nl}^{\pm}\sim(\ell_{\lambda}^{\pm}/\lambda)(\delta\bm{z}^{\pm}_{\lambda}/V_{a})$.222We
assumed here, implicitly, that the cascade is local in $k_{\perp}$. However,
it’s important to note that this condition may not hold in cases of strongly
imbalanced turbulence (Schekochihin, 2022). This point is further discussed in
Section 6.

In the limit of weak, turbulence, $\chi\ll 1$, energy transfer to smaller
scales occurs fractionally upon each collision, with weak nonlinear effects
accumulating over timescales significantly longer than the wave period. Under
the assumption that collisions accumulate in a manner akin to a standard
random walk, approximately $N\sim 1/\chi^{2}$ collisions are necessary to
fully cascade the energy at a given scale. Dimensional analysis (Ng &
Bhattacharjee, 1996) and the formal application of weak perturbative theory
(Galtier et al., 2000) both predict an inertial range spectrum scaling as
$E\propto k_{\perp}^{-2}$.

An initial state of globally weak turbulence is often unstable, and the
intrinsic anisotropy of the energy cascade inevitably develops smaller
perpendicular scales that are strongly turbulent with nonlinear effects
present at the leading order. The case of strong turbulence was addressed by
Higdon (1984) and Goldreich & Sridhar (1995), hereafter GS95. GS95 expanded
upon the idea that fluctuations in any two planes perpendicular to the mean
field can remain correlated only if an Alfvén wave can propagate between them
in less time than their perpendicular decorrelation time (see also, Maron &
Goldreich, 2001; Schekochihin et al., 2009a). This leads to the implication
that the dynamics of the inertial range in incompressible MHD turbulence are
governed by wavevector modes in a “critical balance” (CB) state, i.e.,
characterized by a near-equal balance between the two dynamically important
timescales, essentially achieving $\chi\sim 1$. Under the assumption of
negligible cross helicity, which suggests identical statistical properties for
counterpropagating wavepackets, this implies that the parallel and
perpendicular wavevectors follow the scaling law $k_{||}\propto
k_{\perp}^{2/3}$. Consequently, the one-dimensional power spectra for total
energy scale as $E(k_{\perp})\propto k_{\perp}^{-5/3}$ and $E(k_{||})\propto
k_{||}^{-2}$, across and along the local magnetic field, $\bm{B}_{\ell}$,
respectively.

The GS95 model is directly applicable to the Reduced MHD approximation (RMHD),
where the background magnetic field $\bm{B}_{0}$ is significantly stronger
than the fluctuating amplitudes (Kadomtsev & Pogutse, 1967; Strauss, 1976;
Oughton et al., 2017). The latter are restricted to a plane orthogonal to
$\bm{B}_{0}$. Despite $\delta\bm{b}_{\perp}\ll\bm{B}_{0}$, nonlinear effects
are retained at the leading order. This is achieved by excluding all high-
frequency fluctuations $\tau_{a}\leq\tau_{nl}$. Consequently, fluctuations
within the RMHD approximation inherently satisfy the condition $\chi\geq 1$
(Oughton & Matthaeus, 2020). In this case, nonlinear interactions along the
mean field may be completely neglected, allowing the nonlinear evolution to
adhere to the 2D incompressible MHD equations in planes orthogonal to
$\bm{B}_{0}$333Note, however, that the RMHD approximation encompasses
essential elements of the physics of three-dimensional incompressible MHD
(Dmitruk et al., 2005; Oughton et al., 2017).

While the core principle of the GS95 model, namely CB, was shown to be
consistent with numerical simulations of homogeneous, (in)compressible MHD
turbulence, the inertial range scaling perpendicular to $\bm{B}_{\ell}$ was
observed to be closer to -3/2 (Maron & Goldreich, 2001; Müller et al., 2003;
Müller & Grappin, 2005)444Recent higher resolution simulations further support
this finding (Perez et al., 2012; Verdini & Grappin, 2015; Mallet et al.,
2016; Dong et al., 2022; Shi et al., 2023). Additionally, numerical
simulations revealed a tendency for magnetic and velocity fluctuations in the
field-perpendicular plane to align with each other within a small, scale-
dependent angle (Beresnyak & Lazarian, 2006).

To reconcile the noted discrepancy, Boldyrev (2005, 2006), henceforth B06,
proposed a phenomenological model linking the emergence of local imbalance
(Dobrowolny et al., 1980; Matthaeus et al., 2008) with the scale-dependent
dynamic alignment (SDDA) in the polarizations of $\delta\bm{b}_{\perp}$ and
$\delta\bm{v}_{\perp}$, towards smaller scales, $\theta_{\perp}^{ub}\sim\delta
b/v_{A}\propto\lambda^{1/4}$. Drawing on geometrical considerations, B06
suggests that the observed increase in alignment at smaller scales is linked
to both a depletion of nonlinearities and simultaneous development of local
anisotropy in the plane perpendicular to $\bm{B}_{\ell}$. In this framework
turbulent eddies are identified as 3D-anisotropic structures, characterized by
$\ell_{||}\gg\xi\gg\lambda$, where, $\lambda/\xi\sim\sin\theta_{\perp}^{ub}$.
Here, $\xi$ represents the coherence length in the direction of
$\delta\bm{b}$. B06 predicts three distinct scaling laws in the inertial
range: $E(k_{\xi})\propto k_{\xi}^{-5/3}$, $E(k_{\lambda})\propto
k_{\lambda}^{-3/2}$, and $E(k_{\ell_{||}})\propto k_{\ell_{||}}^{-2}$. The B06
model has been subject to criticism by Beresnyak (2011), on the basis that it
conflicts with the rescaling symmetry intrinsic to the RMHD equations.

The concept of SDDA received further refinement in the work of Chandran et al.
(2015), hereinafter referred to as CSM15. In their interpretation, alignment
emerges as an intermittency effect, resulting from the mutual shearing of
Elsässer fields during the imbalanced collisions ($\delta z^{\pm}\gg\delta
z^{\mp}$) of counterpropagating wave packets. Furthermore, Mallet &
Schekochihin (2017), subsequently referred to as MS17, formulated a
statistical model for three-dimensional 3D anisotropic, intermittent Alfvénic
turbulence. Both models integrate the principles of CB and SDDA, preserving
the scale invariance characteristic of RMHD. With the incorporation of SDDA,
the nonlinear timescale is defined as $\tau_{nl}^{\pm}\sim\lambda/(\delta
z^{\mp}\sin\theta^{z})$, wherein $\theta^{z}$ denotes the angle between the
fluctuations of the Elsässer variables in the plane perpendicular to the
magnetic field, $\delta z^{\pm}_{\perp}$. MS17 extends predictions for 3D
anisotropic higher-order scaling that are consistent with the second-order
moment scaling proposed by B06.

The aforementioned theoretical models focus on the dynamics of small amplitude
(toroidal) Alfvén modes, neglecting potential couplings with compressive
fluctuations (e.g., Cho & Lazarian, 2003; Chandran, 2005, 2018), within the
context of globally balanced turbulence. In contrast, solar wind turbulence
typically displays a predominance in the flux of outwardly propagating
fluctuations over inwardly directed ones (Roberts et al., 1987; D’Amicis et
al., 2021). In addition, the predominant fluctuations in the solar wind are
consistent with large amplitude ($\delta\bm{b}\sim\bm{B}_{0}$) Alfvén waves.
These waves are not purely transverse to $\bm{B}_{0}$ but exhibit spherical
polarization, meaning the magnetic field vector traces a sphere of constant
radius $|\bm{B}_{0}+\delta\bm{b}|=const.$ (Barnes & Hollweg, 1974; Bruno et
al., 2001; Matteini et al., 2014). Furthermore, a non-negligible fraction of
compressive fluctuations is observed (Yao et al., 2011; Howes et al., 2012;
Klein et al., 2012). Finally, the dynamics of the solar wind are influenced
not just by the mean field direction, but also by the radial axis along which
the solar wind expands (Völk & Aplers, 1973). This aspect has been illuminated
by numerical simulations using the Expanding Box Model (EBM) (Grappin et al.,
1993; Grappin & Velli, 1996), which demonstrate that expansion preferentially
reduces the radial component of the magnetic field across all scales,
confining fluctuations to a plane orthogonal to the radial direction (Dong et
al., 2014). Given these considerations, the extent to which phenomenological
models of homogeneous, Alfvénic turbulence can accurately capture the unique
characteristics of the solar wind remains an active area of debate (Bowen et
al., 2021).

In this investigation, we endeavor to evaluate the consistency of the
homogeneous models of balanced MHD turbulence proposed by CSM15 and MS17
against in-situ observations sampled during the first perihelion (E1) of the
Parker Solar Probe mission, (PSP, Fox et al., 2016). Our primary objective is
to conduct a rigorous comparison between the predicted scalings of higher-
order moments in these models and the corresponding empirical observations,
with a specific focus on determining the presence and measurable impact of
model-specific elements, such as SDDA and CB, on the observed characteristics.
Through this analysis, our aim is to not only provide deeper insights into the
higher-order statistics of magnetic turbulence but also to establish a robust
benchmark for the testing and refinement of theoretical models addressing
imbalanced MHD turbulence in the inhomogeneous solar wind.

The structure of the remainder of the paper is as follows: Section 2 provides
theoretical background on intermittency and summarizes the core principles of
the CSM15 and MS17 models. Section 3 elaborates on the methodologies utilized
in this analysis. Details regarding data selection and processing are outlined
in Section 4. The study’s findings are presented in Section 5. A comprehensive
comparison with prior theoretical, observational, and numerical studies, which
contextualizes our results, is provided in Section 6. The paper concludes with
a summary of the main findings in Section 7.

##  2 Theoretical background: Intermittency

Kolmogorov (1941) (hereafter K41), postulated that the cascade process is
intricate enough for eddies to lose all memory of their past and that their
properties after each cascade step can be explained by random distributions.
An inertial range arises where eddies are too large for viscosity to play a
significant role and too small to retain any influence from large-scale
inhomogeneities. A universal probability density function (PDF), $P(\delta
u,\ell)$, emerges over this range. The velocity difference between non-
proximate points can be expressed as the sum of velocity differences over
subintervals, justifying a plausible Gaussian distribution assumption based on
the central limit theorem. The validity of this assumption hinges on two
conditions: mutual independence among summands and comparable finite variances
in the subinterval probability distributions (Feller, 1968). Considering a
global energy transfer rate that is independent of scale, $\epsilon$, and
utilizing the statistical moments of the PDFs called structure functions
($SF^{n}$):

$SF_{i}^{n}(\ell)=\int_{-\infty}^{\infty}{(\delta{u}_{i})^{n}P(\delta{u}_{i},\ell)d(\delta{u}_{i})},$
(2)

where, $\delta{u}_{i}$ denotes longitudinal velocity increments, $\delta
u_{i}=V_{i}(\bm{r}+\bm{\ell})-V_{i}(\bm{r})$, Kolmogorov’s similarity
hypothesis yields:

$SF^{n}_{i}(\ell)\sim\epsilon^{n/3}\ell^{\zeta^{K41}(n)},$ (3)

where $\zeta^{K41}(n)=n/3$, implying global scale invariance (self-similarity)
of the fluctuations.

In the case of MHD turbulence, IK introduced a phenomenological model based on
the coupling of small-scale velocity and magnetic fluctuations by the large-
scale background magnetic field, $\bm{B}_{0}$. Compared to K41, the spectral
transfer is reduced by a factor of $\tau_{A}/\tau_{nl}$ owing to the limited
strength of the interaction exhibited by Alfvén waves. Also in this case the
scaling exponents of $SF^{n}$ are expected to scale linearly with order $n$,
i.e., $\zeta^{IK}(n)=n/4$.

However, it was soon realized experimentally that PDFs of fluctuations in
hydrodynamic (HD) and MHD turbulence tend to display increasingly non-Gaussian
behavior at progressively smaller scales (Batchelor et al., 1949; Burlaga,
1991; Sorriso-Valvo et al., 1999). Additionally, it has been noted that the
spatial inhomogeneity of energy dissipation is expected to alter the scaling
exponents of field increments with respect to length scales $\ell$ (Oboukhov,
1962; Kolmogorov, 1962). Specifically, if the energy transfer rate from larger
to smaller eddies statistically varies with $\ell$, then the dissipation rate
$\epsilon$ in Equation 3 should be substituted with $\epsilon_{\ell}$ (Landau
& Lifshitz, 1959):

$SF^{n}_{i}(\ell)\sim\langle\epsilon_{\ell}^{n/3}\rangle\ell^{n/3}.$ (4)

Expressing $\epsilon_{\ell}^{n/3}$ through a scaling relation with $\ell$, we
find $\langle\epsilon_{\ell}^{n/3}\rangle\sim\ell^{\tau_{n}/3}$, leading to

$SF^{n}_{i}(\ell)\sim\ell^{\zeta_{n}},$ (5)

where, $\zeta_{n}=n/3+\tau_{n}/3$ is generally a nonlinear function of n. When
analyzing the scaling exponents $\zeta_{n}$ of the structure functions,
deviations from linear scaling, indicate a violation of global scale
invariance, which implies a process characterized by “multifractal” statistics
and intermittency, i.e., the concentration of the energy into smaller volumes
of space at smaller scales (Frisch, 1995). Several authors have put forward
theoretical arguments to account for the intermittency effect (Gurland &
Tripathi, 1971; Frisch et al., 1978; She & Leveque, 1994; Grauer et al., 1994;
Politano & Pouquet, 1995; Ruzmaikin et al., 1995; Horbury & Balogh, 1997a;
Müller & Biskamp, 2000; Boldyrev, 2002).

Chandran, Schekochihin, and Mallet Chandran et al. (2015) (henceforth CSM15)
introduced a model founded on collisions of Alfvénic wavepackets making a
series of plausible assumptions regarding the dynamics and statistics of RMHD
turbulence. In formulating an analytical model, they proposed two primary
types of non-linear interactions. The first involves occasional “balanced”
interactions, $\delta\bm{z}^{-}\approx\delta\bm{z}^{+}$, where the amplitude
of the wavepackets is reduced by a factor of $0\leq\beta\leq 1$, yet their
length scale remains unchanged. Consequently, the amplitude of a fluctuation
at scale $\lambda$ resulting from such balanced collisions can be expressed as
follows:

$\delta\bm{z}^{\pm}_{\lambda}=\delta\bar{z}\beta^{q}.$ (6)

Here, $\delta\bar{z}$ is the initial fluctuation amplitude at the injection
scale $L$, and $q$ is the number of balanced collisions expected as the
fluctuation evolves from scale $L$ to $\lambda$. The value of $q$ is presumed
to follow a Poisson distribution:

$P(q)=\frac{e^{-\mu}\mu^{q}}{q!},$ (7)

where $\mu$ is the scale-dependent mean value of $q$. The “typical”
fluctuation amplitude that best characterizes the bulk of the volume is

$\delta\bm{z}^{\ast}_{\lambda}=\delta\bar{z}\beta^{\mu}.$ (8)

On the contrary, for “imbalanced”, $\delta\bm{z}^{\pm}\gg\delta\bm{z}^{\mp}$,
collisions, the amplitudes of the fluctuations remain constant while the sub-
dominant field is sheared into alignment and its perpendicular scale $\lambda$
reduces. Assuming that the most intense coherent structures in MHD turbulence
–specifically, those with $q=0$– are 3D anisotropic current sheets with a
volume filling factor $f_{cs}\propto\lambda$, this leads to a straightforward
relation for the scaling exponents:

$\zeta^{CSM15}(n)=1-\beta^{n}.$ (9)

By focusing solely on the scenario of zero cross helicity, CSM15 deduce a
value for $\beta\approx 0.691$. This implies that $\zeta_{n}\rightarrow 1$ as
$n\rightarrow\infty$.

Furthermore, the model offers predictions for the alignment angles between the
perpendicular components of velocity-magnetic and/or Elsässer field
fluctuations defined in various ways. Conventionally, the alignment is
estimated as

$sin(\theta^{ub}_{\perp})=\left\langle\frac{|\delta\bm{b}_{\perp}\times\delta\bm{v}_{\perp}|}{|\delta\bm{b}_{\perp}||\delta\bm{v}_{\perp}|}\right\rangle.$
(10)

A similar expression can be formulated for the angle between Elsässer field
fluctuations, $\delta\bm{z}^{\pm}_{\perp}$, which relates to the residual of
magnetic over kinetic energy, usually studied through the lens of the
normalized residual energy :

$\sigma_{r}=\frac{\langle\delta\bm{v}_{\perp}^{2}\rangle-\langle\delta\bm{b}_{\perp}^{2}\rangle}{\langle\delta\bm{v}_{\perp}^{2}\rangle+\langle\delta\bm{b}_{\perp}^{2}\rangle}.$
(11)

Similarly, the Elsässer imbalance is assessed by comparing the relative energy
in inwardly and outwardly propagating Alfvén waves (Velli et al., 1991; Velli,
1993):

$\sigma_{c}=\frac{2\langle\delta\bm{v}_{\perp}\cdot\delta\bm{b}_{\perp}\rangle}{\langle\delta\bm{v}_{\perp}^{2}\rangle+\langle\delta\bm{b}_{\perp}^{2}\rangle}$
(12)

Observations from numerical simulations showed little to no scaling of
alignment angles as defined in Equation 10, (Beresnyak & Lazarian, 2009).
However, an alternate definition of these angles, achieved by separately
averaging the numerator and denominator –referred to as polarization
intermittency– was proposed by (Beresnyak & Lazarian, 2006):

$sin(\tilde{\theta}^{ub}_{\perp})=\frac{\langle|\delta\bm{b}_{\perp}\times\delta\bm{v}_{\perp}|\rangle}{\langle|\delta\bm{b}_{\perp}||\delta\bm{v}_{\perp}|\rangle}.$
(13)

In a similar manner, we can define the angle between the two Elsässer fields
as $\tilde{\theta}^{z}_{\perp}$. The model by CSM15 predicts
$\tilde{\theta}^{ub}_{\perp}\propto\lambda^{0.21}$ and
$\tilde{\theta}^{z}_{\perp}\propto\lambda^{0.10}$.

Mallet & Schekochihin (2017) (hereafter, MS17) formulated a statistical model
of RMHD turbulence grounded in three principal concepts: critical balance,
dynamic alignment, and intermittency. To substantiate their model, they put
forth four conjectures based on physical reasoning: (a) the fluctuation
amplitudes adhere to an anisotropic log-Poisson distribution (Chandran et al.,
2015; Zhdankin et al., 2016); (b) the structures at small scales are 3D
anisotropic with a sheet-like morphology (Boldyrev, 2006; Howes, 2015); (c)
the scale-independence of the critical balance parameter, inclusive of dynamic
alignment, within the inertial range (Mallet et al., 2015); (d) a consistent
energy flux across parallel scales in the inertial range (Beresnyak, 2015).
The model offers predictions for scaling in the perpendicular, parallel, and
fluctuation directions:

For the perpendicular direction:

$\zeta^{MS17}_{\lambda}(n)=1-\beta^{n},$ (14)

for the parallel direction:

$\zeta^{MS17}_{\ell_{||}}(n)=2(1-\beta^{n}),$ (15)

and for the fluctuation direction:

$\zeta^{MS17}_{\xi}(n)=\frac{n(1-\beta^{n})}{n/2+1-\beta^{n}},$ (16)

where $\beta=1/\sqrt{2}$. A more practical discussion of the model is
presented in Schekochihin (2022).

##  3 Data Analysis

To investigate and quantify the three-dimensional anisotropy of higher-order
magnetic field moments, we employ a methodology proposed by Wang et al.
(2022), which builds upon and extends the framework established in Chen et al.
(2012).

Adhering to the approach outlined in Chen et al. (2012), we establish a
locally-defined, scale-dependent Cartesian coordinate system, represented as
($\hat{\xi}$, $\hat{\lambda}$, $\hat{\ell_{||}}$). In this coordinate system
the “parallel” direction, $\hat{\ell_{||}}$, is aligned with the local
magnetic field, $\bm{B_{\ell}}$, defined by:

$\bm{B_{\ell}}=[\bm{B}(\bm{r}+\bm{\ell})+\bm{B}(\bm{r})]/2,$ (17)

where $\bm{\ell}=\tau(\bm{V}_{\ell}-\bm{V}_{\text{sc}})$ represents the
displacement vector, where $\bm{V}_{\ell}$ and $\bm{V}_{\text{sc}}$ are the
local scale-dependent velocity field and the spacecraft velocity,
respectively, in the RTN coordinate system (Franz & Harper, 2002). Magnetic
field increments are calculated using:

$\delta\bm{b}=\bm{B}(\bm{r}+\bm{\ell})-\bm{B}(\bm{r}).$ (18)

The amplitude of the field increment is denoted as $\delta b=|\delta\bm{b}|$.
The local “displacement” direction, $\hat{\xi}$, aligns with the unit vector
of the perpendicular component of the field increment, with
$\delta\bm{b}_{\perp}$ defined as:

$\delta\bm{b}_{\perp}=\bm{B_{\ell}}\times(\delta\bm{b}\times\bm{B_{\ell}}).$
(19)

Lastly, the “perpendicular” direction, $\hat{\lambda}$, is orthogonal to both
$\hat{\xi}$ and $\hat{\ell_{||}}$,
$\hat{\lambda}=\hat{\ell_{||}}\times\hat{\xi}$. The Cartesian system can be
converted into a spherical polar coordinate system ($\ell$, $\theta_{B}$,
$\phi_{\delta\bm{B}_{\perp}}$), where $\theta_{B}$ is the angle between
$\bm{B}_{\ell}$ and $\bm{\ell}$, and $\phi_{\delta\bm{B}_{\perp}}$ is the
angle between $\hat{\xi}$ and the projection of $\bm{\ell}$ onto the plane
orthogonal to $\bm{B}_{\ell}$.

In our analysis, we utilize the 5-point (5-point) increment method, which
represents a significant advancement over the conventional 2-point (2-point)
method, especially for examining turbulence statistics in sub-ion regimes
(Cerri et al., 2019). A critical benefit of the 5-point method is its reduced
susceptibility to large-scale spectral leakage (Cho, 2019), rendering it more
suitable and effective for our ensuing analysis.

For calculating 5-point structure functions, denoted as $SF^{n}_{5}$, a
modified definition of $\delta\bm{b}$ is required

$\displaystyle\delta\bm{b}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{35}}\left[\bm{B}(\bm{r}-2\bm{\ell})-4\bm{B}(\bm{r}-\bm{\ell})\right.$
(20)
$\displaystyle\left.+6\bm{B}(\bm{r})-4\bm{B}(\bm{r}+\bm{\ell})+\bm{B}(\bm{r}+2\bm{\ell})\right].$

Moreover, the local scale-dependent value $\bm{\psi}_{\ell}$ of a field
$\bm{\psi}$ can be computed as a weighted average using a five-point stencil:

$\displaystyle\bm{\psi}_{\ell}$ $\displaystyle=$
$\displaystyle\frac{1}{16}\left[\bm{\psi}(\bm{r}-2\ell)+4\bm{\psi}(\bm{r}-\ell)\right.$
(21)
$\displaystyle\left.+6\bm{\psi}(\bm{r})+4\bm{\psi}(\bm{r}+\ell)+\bm{\psi}(\bm{r}+2\ell)\right].$

For example, the local scale-dependent magnetic and velocity fields are
represented by $\bm{B}_{\ell}$ and $\bm{V}_{\ell}$, respectively.

The $\mathit{n}$th-order, structure functions conditioned on the pair of
angles $\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}}$, are defined as:

$SF^{n}(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})=\langle(\delta
B)^{n}|~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}},~{}\ell\rangle,$ (22)

The conditional average in Equation 22 was calculated over the angle bin
$\omega(i-1)^{\circ}\leq\theta_{B}\leq\omega i^{\circ}$,
$\omega(j-1)^{\circ}\leq\phi_{\delta\bm{B}_{\perp}}\leq\omega j^{\circ}$,
where $i=1,...,9$ and $j=1,...,9$. In the following, $\omega$ takes the value
of $\omega=5$ for estimating lower-order moments and $\omega=10$ for higher-
order moments.

In the following, we focus mainly on three special cases, defining the
components in

$i=1,\quad j=1:\quad SF(\ell_{||})^{n},\quad\quad\text{``parallel''},\quad$
(23) $i=9,\quad j=1:\quad SF(\xi)^{n},\quad\text{``displacement''},$ (24)
$i=9,\quad j=9:\quad SF(\lambda)^{n},\quad\text{``perpendicular''},$ (25)

directions, where
$\ell_{||}=\bm{\ell}\cdot\bm{\hat{\ell}_{||}},~{}\lambda=\bm{\ell}\cdot\bm{\hat{\lambda}}$,
and $\xi=\bm{\ell}\cdot\bm{\hat{\xi}}$.

To estimate a component structure function for the entire dataset, we adopt
the methodology outlined in Verdini et al. (2018). For each selected interval,
$j$, within the dataset, we compute the structure functions
$SF^{n}_{j}(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})$ for the three
orthogonal components, as defined by Equations 23 to 25. Considering the
substantial variation in the root mean square (rms) of fluctuations between
intervals, normalization is a critical step prior to averaging these
intervals. The normalization involves selecting a specific scale,
$\ell_{\ast}$, and normalizing each $SF^{n}_{j}$ by the energy of fluctuations
at that scale. We determine an appropriate $\ell_{\ast}$ by estimating the
trace structure function $S^{n}_{j}(\ell)$ for each interval and identifying a
scale range where power-law behavior is consistent across all
$S^{n}_{j}(\ell)$. The fluctuation energy at scale $\ell_{\ast}$ is given by
the value of the trace structure function $S_{i}(\ell_{\ast})$.

The normalized weighted average structure function for a given magnetic field
component is then calculated as follows:

$\tilde{SF}^{n}(\ell,\theta_{B},\phi_{\delta\bm{B}_{\perp}})=\sum_{j}\frac{SF^{n}_{j}(\ell,\theta_{B},\phi_{\delta\bm{B}_{\perp}})}{S^{n}_{j}(\ell_{\ast})}W_{j},$
(26)

Figure 1: The local 2-point structure functions (circles) and 5-point
structure functions (squares) averaged for the five most highly Alfvénic
intervals within our dataset. The structure functions are displayed for the
parallel, perpendicular, and displacement directions, indicated by red, black,
and blue colors, respectively. Reference lines representing scalings of 1/2,
2/3, and 1 are included for comparison. Figure 2: Overview of E1: (a) Magnetic
field timeseries, radial component, $B_{R}$ (blue) and magnitude, $|B|$
(black); (b) Solar wind speed, $V_{sw}$ (black, left axis) and proton
temperature, $T_{p}$ (blue, right axis); (c) Normalized cross helicity,
$\sigma_{c}$ (black) and normalized residual energy, $\sigma_{r}$ (red); (d)
Plasma $\beta$ (black, left axis) and Alfvénic Mach number,
$M_{a}=V_{R}/|V_{a}|$ (blue, right axis); (e) Variance anisotropy,
$E=(b^{2}_{T}+b^{2}_{N})/b^{2}_{R}$, where $b$ represents the rms amplitude of
fluctuations (black, left axis) and angle between the magnetic field and
velocity flow, $\Theta_{VB}$ (blue); (f) Sampling angle, $\Theta_{R}$, defined
as the angle between $\hat{R}$ and $\bm{V}_{sc}-\bm{V}_{sw}$ (black, left
axis), and radial distance from the Sun, $R$ (red, right axis). Additionally,
three intervals denoted as $I_{j},~{}j=~{}1,~{}2,~{}3$ and marked with black,
pink, and cyan shadings on the main figure. The corresponding 5-point
structure functions of the parallel, perpendicular, and displacement
directions, denoted by red, transparent black, and blue colors, respectively,
are shown in panels (g)-(i). In addition, the trace structure function is
shown in black circles. Reference lines representing scalings of 1/2, 2/3, and
1 are included for comparison.

where
$W_{j}=n_{j}(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})/n(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})$
represents the weighting factor, with
$n_{j}(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})$ being the total
number of measurements within each bin for the interval under consideration,
divided by the count in each bin for the whole data set, $n=\sum_{j}n_{j}$.

Our subsequent analysis, as discussed in Section 4, relies on these
conditionally defined structure functions, utilizing data from the first
perihelion of the Parker Solar Probe (PSP). Unless specified otherwise, the
results in the following analysis are derived from estimating 5-point
structure functions.555The algorithm detailed in this section, along with a
package for downloading, cleaning, and processing PSP data, is readily
accessible in MHDTurbPy (Sioulas, 2023).

##  4 Data Set

We analyze magnetic field and particle data collected during the first
perihelion of the PSP mission, covering the period from November 1 to November
11, 2018. Specifically, we analyze magnetic field measurements obtained by the
FIELDS instrument (Bale et al., 2016). In particular, we make use of the SCaM
data product, which combines measurements from fluxgate and search-coil
magnetometers (SCM) by using frequency-dependent merging coefficients. This
approach allowed us to observe the magnetic field over a frequency range
ranging from direct current (DC) to 1 MHz while achieving optimal signal-to-
noise ratio (Bowen et al., 2020). The FIELDS magnetometer suite is susceptible
to narrow-band coherent noise stemming from the spacecraft reaction wheels,
including their rotation frequencies, as well as harmonic and beat
frequencies. To address potential contamination of magnetic field measurements
at ion-scales by the reaction wheels, for each interval, we implement a
procedure involving the identification and elimination of reaction wheel noise
using the method detailed in Shankarappa et al. (2023).

Moreover, we incorporated data from the Solar Probe Cup (SPC) instrument,
which is part of the Solar Wind Electrons Alphas and Protons (SWEAP) suite
(Kasper et al., 2016), to estimate bulk plasma properties. We also utilized
Quasi Thermal Noise (QTN) electron density measurements (Moncuquet et al.,
2020; Pulupa et al., 2017). To enhance statistical robustness and expand the
sample size, the data were segmented into 12-hour intervals. These intervals
were designed to overlap by 6 hours to maximize data utilization. We then
conditioned the intervals based on $\sigma_{c}$, selecting only those with an
average value of $\sigma_{c}(\ell_{\ast})\geq 0.75$, where
$\ell_{\ast}=10^{4}d_{i}$. This methodology yielded a total of 82 intervals,
sampled at distances ranging from 0.166 to 0.244 au.

Table 1:  The median values for the spectral indices of the trace, parallel,
perpendicular, and displacement components of the magnetic field in the ranges
$R_{1}$ and $R_{2}$. These indices are derived from the corresponding scaling
indices of $SF_{5,~{}i}^{2}$, utilizing the relationship
$\alpha_{i}=-1-\beta_{i}$ (Monin & Jaglom, 1987). The error values provided
represent the standard deviation of the mean.

| $\alpha$ | $\alpha_{\ell_{||}}$ | $\alpha_{\lambda}$ | $\alpha_{\xi}$
---|---|---|---|---
$R_{1}$ | $-1.79\pm 0.06$ | $-1.97\pm 0.05$ | $-1.64\pm 0.04$ | $-1.94\pm 0.06$
$R_{2}$ | $-1.53\pm 0.02$ | $-1.66\pm 0.05$ | $-1.49\pm 0.03$ | $-1.56\pm 0.08$

Figure 3: Wavevector anisotropy relationships: (a)
$\ell_{||}=\ell_{||}(\lambda)$ (b) $\xi=\xi(\lambda)$, determined by equating
pairs of $SF^{2}_{5}(\lambda)$ with $SF^{2}_{5}(\ell_{||})$, and
$SF^{2}_{5}(\lambda)$ with $SF^{2}_{5}(\xi)$, respectively. The gray lines
show the aspect ratios, $\ell_{||}/\lambda$ and $\xi/\lambda$ plotted against
$\lambda$ in panels (a) and (b), respectively. For context, reference lines
indicating scalings of 2/3 and 1 are also included. Panel (c) presents power
anisotropy, with $E_{\lambda}/E_{\ell_{||}}$ depicted in red and
$E_{\lambda}/E_{\xi}$ in blue.

##  5 Results

###  5.1 A Comparison of $SF_{5}^{n}$ and $SF_{2}^{n}$

This section is dedicated to a comparative analysis aimed at substantiating
our preference for the 5-point structure function method over the
traditionally favored 2-point approach. To this end, we calculated second-
order structure functions for the parallel, perpendicular, and displacement
components according to Equations 23 to 25, setting $\omega=5$. This involves
estimating $\bm{B}_{\ell}$ and $\delta\bm{b}$ as per Equations 17 to 19 for
$SF_{2}^{2}$, and Equations 21 to 22 for $SF_{5}^{2}$.

The local 2-point structure functions (circles) and 5-point structure
functions (squares) averaged for the five most highly Alfvénic intervals
within our dataset are illustrated in Figure 1. To highlight specific scale
ranges, pink and gray shadings are employed for intervals $8-100d_{i}$
(labeled $R_{1}$) and $200-6000d_{i}$ (labeled $R_{2}$), respectively. At
large scales, comparable results are obtained from both the 5-point and
2-point methods. However, a marked divergence is observed in the $R_{1}$
range. The 5-point method reveals steep scaling for parallel and displacement
components, with indices $\beta_{\ell_{||}}\approx\beta_{\xi}\approx 1$,
aligning with wavelet-derived parallel component scaling of the same dataset
reported in Sioulas et al. (2023). In contrast, the 2-point method produces a
notably flatter slope.

The distinction between the two methods becomes stark at kinetic scales,
highlighting the 2-point method’s limitations in detecting steep scalings.
This shortfall is further evident when comparing $\bm{B}$-trace wavelet
structure functions using both the $SF_{2}^{2}$ and $SF_{5}^{2}$ methods; only
the 5-point approach produces scalings that align with wavelet analyses across
scales from injection to kinetic (details not shown here). This inconsistency
underscores potential inaccuracies when employing the 2-point method in
scenarios characterized by steep scaling, as also emphasized by Cerri et al.
(2019)666A cautionary note is warranted: when comparing $SF_{5}^{2}$ with
wavelet-derived trace structure functions for intervals observed in later PSP
encounters —characterized by shallower than $1/f$ energy injection scale power
spectra (Huang et al., 2023; Davis et al., 2023)— $SF_{5}^{2}$ fails to
replicate the wavelet scalings. This indicates that $SF_{5}^{2}$, similar to
its $SF_{2}^{2}$ counterpart, is ineffective under conditions with scalings
shallower than $1/f$.. Therefore, for our further analyses, we have chosen to
rely on the 5-point method.

Across all examined intervals, the results remained qualitatively consistent.
Three such intervals are illustrated in panels (g)-(h) of Figure 2. In
addition to the component structure functions, these panels also illustrate
the trace structure function estimated for the respective intervals. It can be
observed that the perpendicular component can significantly diverge from the
trace, while the latter typically shows remarkable overlap with the
displacement component, with the difference becoming more pronounced towards
smaller scales.

Table 1 presents the median $SF_{5}^{2}$ scalings estimated across the $R_{1}$
and $R_{2}$ ranges for the entire dataset.

###  5.2 Power & wavevector anisotropy

To examine the scale-dependent three-dimensional anisotropy in our dataset, we
calculated $SF_{5}^{2}(\ell)$, following the methodology outlined in Equations
23 to 25, setting $\omega=5$. Our analysis commenced by identifying the
anisotropic relationships for individual intervals, then proceeded to compute
a scale-dependent median for the entire dataset, utilizing 150 logarithmically
spaced bins.

Panels (a) and (b) of Figure 3 illustrate the anisotropic relationships
$\ell_{||}(\lambda)$, and $\xi(\lambda)$, derived by equating pairs of
structure functions: $SF^{2}_{5}(\lambda)$ with $SF^{2}_{5}(\ell_{||})$, and
$SF^{2}_{5}(\lambda)$ with $SF^{2}_{5}(\xi)$, respectively. The aspect ratios
$\ell_{||}/\lambda$ and $\xi/\lambda$ are represented by gray lines in their
respective panels. Panel (c) focuses on power anisotropy, illustrating the
ratios $SF^{2}_{5}(\lambda)/SF^{2}_{5}(\ell_{||})$ and
$SF^{2}_{5}(\lambda)/SF^{2}_{5}(\xi)$ in red and blue, respectively. The
median values of the dataset’s anisotropic scalings are summarized in Table 2.

At large scales, our observations reveal a rough equipartition of the
fluctuating energy between $SF^{2}_{5}(\lambda)$ and $SF^{2}_{5}(\ell_{||})$,
with $SF^{2}_{5}(\xi)$ being slightly more energetic. The energy distribution
is reflected in the wavevector anisotropy and aspect ratios, indicating that
eddies tend to be slightly compressed along the fluctuation direction.

Table 2:  Median values of the scaling indices for wavevector anisotropies,
specifically $\ell_{||}\propto\lambda^{w_{\ell_{||}}}$ and
$\xi\propto\lambda^{w_{\xi}}$, and the power anisotropies
$E_{\lambda}/E_{\ell_{||}}\propto\lambda^{p_{\ell_{||}}}$ and
$E_{\lambda}/E_{\xi}\propto\lambda^{p_{\xi}}$. These indices were derived by
applying a power-law fit to the curves obtained from individual intervals over
the scale ranges $R_{1}$ and $R_{2}$. The table presents the median values
along with their associated errors, represented as the standard deviation of
the mean.

| $w_{\ell_{||}}$ | $w_{\xi}$ | $p_{\ell_{||}}$ | $p_{\xi}$
---|---|---|---|---
$R_{1}$ | $0.89\pm 0.06$ | $0.77\pm 0.05$ | $-0.31\pm 0.04$ | $-0.30\pm 0.07$
$R_{2}$ | $0.86\pm 0.08$ | $0.99\pm 0.06$ | $-0.19\pm 0.05$ | $-0.09\pm 0.06$

Figure 4: 3D representation of turbulence eddies obtained by estimating
isosurfaces of constant energy of $SF^{2}_{5}$ at different levels, ranging
from small scales (top left) to large scales (bottom right). The color scheme,
although redundant, indicates the distance from the origin (0,0,0).
Additionally, projections of the object onto each respective plane are
displayed. While the projections share a common colormap to denote the
distance from each plane’s origin, the colormap for the 3D object differs.

Within the $R_{2}$ range, we note that the fluctuating magnetic energy is
distributed almost evenly between the perpendicular and displacement
components. This observation suggests rough axisymmetry of the turbulent
eddies at these scales, i.e., approximate isotropy in the plane perpendicular
to $\bm{B}_{\ell}$. The aspect ratio $\xi/\lambda$ is observed to be close to,
yet marginally greater than, unity, displaying only a slight increase within
$R_{2}$. This trend is further emphasized by the scale-dependent power-
anisotropy depicted in panel (c). In contrast, within this range, eddies
exhibit elongation along $\bm{B}_{\ell}$, as indicated by the ratio
$\ell_{||}/\lambda$, which is greater than 1 and shows a monotonic increase
towards smaller scales. These findings collectively indicate that turbulent
eddies within the $R_{2}$ range predominantly exhibit a field-aligned tube
topology, consistent with the results presented in Vinogradov et al. (2023).
This is visually demonstrated in Figure 4. At scales smaller than
$\lambda\lesssim 100d_{i}$, a noticeable shift from isotropy in the plane
perpendicular to $\bm{B}_{\ell}$ becomes evident. This shift is highlighted by
a gradual increase in the aspect ratio $\xi/\lambda$, indicating a transition
of eddy structures from tube-like to ribbon-like. The evolving eddies exhibit
three-dimensional anisotropy, adhering to the relationship
$\ell_{||}\gg\xi\gg\lambda$. This trend persists into the smaller-scale end of
the $R_{1}$ range, where both power-anisotropy ratios follow a scaling of
approximately -1/3. Additionally, within $R_{1}$, the ratio
$\ell_{||}/\lambda$ continues to rise, albeit at a reduced rate compared to
the $100-600d_{i}$ range, where $\ell_{||}\propto\lambda^{0.72\pm 0.04}$. This
slower rate of increase within $R_{1}$ can be attributed to the steepening of
$SF^{2}_{5}(\ell_{||})$ at scales marginally larger than $R_{1}$.

As we move below the $R_{1}$ range and into the transition region (Sahraoui et
al., 2009; Bowen et al., 2020a), the previously observed trend of increasing
anisotropy ceases. Within the scale range of $2d_{i}\lesssim\lambda\lesssim
8d_{i}$, the eddies start to demonstrate more isotropic characteristics. This
tendency towards isotropy peaks at $\lambda=2d_{i}$, where the aspect ratio
reaches $\lambda:\xi:\ell=1.56:1:1$. Throughout this scale range, the
anisotropic scaling relations —derived from fitting curves to individual
intervals and estimating the median values, presented along with the standard
deviation of the mean— conform to $\ell_{||}\propto\lambda^{2.01\pm 0.06}$ and
$\xi\propto\lambda^{1.25\pm 0.05}$.

At even smaller scales, distinct scaling anisotropies characterize two
separate ranges. Within $1d_{i}\lesssim\lambda\lesssim 2d_{i}$, the ratio
$\ell_{||}/\lambda$ shows an upward trend, following
$\ell_{||}\propto\lambda^{0.67\pm 0.02}$, while the ratio $\xi/\lambda$
remains relatively stable, adhering to $\xi\propto\lambda^{0.98\pm 0.04}$.
These findings stand in contrast to standard kinetic Alfvén wave (KAW)
turbulence models (Howes et al., 2008; Schekochihin et al., 2009b) and deviate
statistically from in-situ observations reported in (Duan et al., 2021; Zhang
et al., 2022). However, they are in agreement with the intermittent KAW model
proposed by (Boldyrev & Perez, 2012) and align with numerical kinetic
simulations by (Cerri et al., 2019), as well as in-situ observations in the
magnetosheath (Wang et al., 2020). At yet smaller scales,
$0.5d_{i}\leq\lambda\leq 1d_{i}$, both $\ell_{||}/\lambda$ and $\xi/\lambda$
ratios exhibit an increase, scaling as $\ell_{||}\propto\lambda^{0.5\pm 0.05}$
and $\xi\propto\lambda^{0.87\pm 0.04}$, respectively.

Figure 4 presents a three-dimensional representation of turbulence eddies,
illustrating isosurfaces of $SF^{5}_{2}$ at various scales. This visualization
was achieved by estimating conditional
$SF^{2}_{5}(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})$ , according to
Equation 22, and utilizing $5^{\circ}$ angular bins. The spherical polar
coordinates $(\ell,~{}\theta_{B},~{}\phi_{\delta\bm{B}_{\perp}})$ obtained
from this process were then converted into Cartesian coordinates
$(\ell_{||},~{}\xi,~{}\lambda)$. Surfaces computed for the first octant were
mirrored across to the other octants, based on the assumption of reflectional
symmetry (Chen et al., 2012). In the 3D visualization, surface colors
represent the distance from the origin, with cooler colors indicating larger
distances. When these surfaces are projected onto different planes, the color
denotes the distance from the origin of each respective plane. It is important
to note that the colormap applied to these planar projections differs from the
one used for the 3D representation. The color coding in these projections
reflects the range of maximum and minimum distances observed across all three
components.

###  5.3 Higher order statistics & Intermittency

We computed five-point structure functions, $SF_{5}^{n}$, for parallel,
perpendicular, and displacement components, as per Equations 23 to 25,
considering orders $n=1,\ldots,10$ with $\omega=10$. Additionally, we
evaluated $\bm{B}$-trace structure functions, $S_{5}^{n}$. Furthering our
analysis, we derived thenormalized weighted average, $\tilde{SF}{5}^{n}$, for
our dataset following Equation 26. After identifying two distinct sub-inertial
ranges displaying clear power-law behavior, we calculated two sets of
normalized structure functions. Within the $R_{1}$ domain, we normalized
$SF_{5}^{n}$ using $S_{5}^{n}(\lambda_{\ast})$, where
$\lambda_{\ast}=50d_{i}$.

Figure 5: The main panels depict the normalized weighted average structure
functions, $\tilde{SF}^{n}_{5}(\ell)$, for (a) parallel, (b) perpendicular,
and (c) displacement components, each adjusted with a vertical offset for
enhanced clarity. Prior to the weighted average estimation, each
$SF^{n}_{5}(\ell)$ was normalized using the value of the trace
$S^{n}_{5}(\lambda{\ast})$, where $\lambda_{\ast}=2\times 10^{3}d_{i}$. These
normalized averages, $\tilde{SF}^{n}_{5}(\ell)$, displayed in the main figure,
were used to estimate $\zeta_{n}$ for $R_{2}$, marked as gray asterisks in the
insets. For $\zeta_{n}$ estimation in $R_{1}$, indicated by red asterisks in
the insets, a similar normalization process was applied at
$\lambda_{\ast}=50d_{i}$. It is important to note that the latter normalized
$\tilde{SF}^{n}_{5}(\ell)$ for $R_{1}$ are not depicted in the figure. Error
bars in the insets represent the uncertainty associated with the power-law
fits. For comparison, the scaling behaviors as predicted by the K41, IK,
CSM15, and MS17 models are also included.

For the $R_{2}$ domain, normalization employed a different scale,
$\lambda_{\ast}=2\times 10^{3}d_{i}$. The latter normalized structure
functions, particularly for the parallel, perpendicular, and displacement
components, are depicted in Figure 5 panels (a) to (c), respectively. To aid
visualization, each n-th order moment was vertically offset by $10^{-n}$.
Fitting each component of $\tilde{SF}_{5}^{n}(\ell)$ to a power law,
$\propto\ell^{\zeta_{n}}$, facilitated the estimation of scaling exponents,
$\zeta_{n}$. The resulting $\zeta_{n}$ are depicted by red asterisks for the
$R_{1}$ domain and gray for $R_{2}$ within the insets of the corresponding
figure panels. For comparison, scaling exponent predictions based on the
theoretical models proposed by CSM15 and MS17 are also included.

In the $R_{2}$ range, the scaling exponent $\zeta_{n}^{\lambda}$ of the
perpendicular component forms a convex function of n, indicative of
multifractal statistics and strong intermittency. This scaling profile closely
aligns with the theoretical predictions by CSM15 and MS17, showing notable
correspondence to the latter model at lower $n$ values. For $SF^{n}_{5}(\xi)$,
the observed scaling exponents, $\zeta_{n}^{\xi}$, slightly deviate towards
shallower gradients compared to the MS17 model. The scaling exponents of the
parallel component $\zeta_{n}^{\ell_{||}}$ exhibit a nonlinear dependence on
$n$, though with less pronounced concavity than the perpendicular components,
and notably diverge from the MS17 model, even at lower $n$ values. A
comprehensive discussion of these findings and their broader implications is
provided in Section 6.

In the small-scale sub-inertial range, $R_{1}$, the scaling exponents for both
the parallel and displacement components display a linear relationship with
$n$. Conversely, the scaling exponent $\zeta_{n}^{\lambda}$ of the
perpendicular component forms a convex function of n, albeit demonstrating a
lesser extent of non-linearity relative to that observed in the $R_{2}$ range.
Overall, $R_{1}$ is characterized by less pronounced intermittency signatures
compared to $R_{2}$, with the observed $\zeta_{n}$ profiles deviating from the
expectations set by established theoretical models.

To further investigate multifractality and deviations from Gaussian statistics
in the magnetic field time series, we consider the Scale-Dependent Kurtosis,
defined as $K(\ell)=SF^{4}(\ell)/[SF^{2}(\ell)]^{2}$ (Frisch, 1995; Bruno et
al., 2003). As a normalized fourth-order moment, $K(\ell)$ is sensitive to
extreme values of increments, allowing us to detect the tendency of PDFs in
intermittency-affected time series to exhibit increasingly flared-out tails at
smaller scales. In simpler terms, it quantifies how the “tailedness” of the
distribution of increments in a turbulent field changes across various scales.

We employ both 2-point ($K_{2}(\ell)$) and 5-point ($K_{5}(\ell)$) methods to
study the fractal properties of magnetic field time series. The resulting
curves for the parallel, perpendicular, and displacement components of the
magnetic field are illustrated in panels (a) to (c) of Figure 6. At scales
$\lambda\geq 100d_{i}$, an increase in $K(\ell)$ is observed for all
components towards smaller scales, indicating a progressive deviation from
Gaussianity in the underlying PDFs of increments, a hallmark of multifractal
statistics (Sorriso-Valvo et al., 1999).

The limitations of the 2-point method are particularly evident in the $R_{1}$
range, where it notably diverges from the 5-point approach. In the case of the
perpendicular component, $K_{2}(\lambda)$ appears to plateau at
$\lambda<20d_{i}$, consistent with the findings of Chhiber et al. (2021).
Conversely, $K_{5}(\ell)$ maintains an increasing trend in $R_{1}$, albeit
with a less steep slope compared to $R_{2}$. Within the $R_{1}$ range, both
the parallel, $K_{5}(\ell_{||})$, and displacement, $K_{5}(\xi)$, components
exhibit super-Gaussian but monofractal behavior, consistent with the linear
$\zeta_{n}^{\ell_{||}}$ and $\zeta_{n}^{\xi}$ profiles illustrated in Figure
5.

At kinetic scales, both $K_{5}(\lambda)$ and $K_{5}(\xi)$ exhibit an
increasing trend towards smaller scales, while the behavior of the parallel
component remains less distinct. This trend in the perpendicular component
diverges from the monofractal statistics observed at sub-ion scales in several
observational studies using $K_{2}$ (Wu et al., 2013; Chen et al., 2014;
Chhiber et al., 2020). However, our findings align with hybrid and fully
kinetic simulations by Cerri et al. (2019), where $K_{5}(\lambda)$ is
demonstrated to increase above Gaussian values throughout the sub-ion scale
range. Furthermore, qualitatively consistent trends were identified by
Alexandrova et al. (2008), who utilized wavelet-derived kurtosis to observe a
gradual increase in this measure at kinetic scales.

###  5.4  Scale-Dependent Dynamic Alignment & Critical Balance

We begin by examining the scale-dependent behavior of the alignment angle
between the perpendicular components of the increments,
$\delta\bm{b}_{\perp}-\delta\bm{u}_{\perp}$ and
$\delta\bm{z}^{+}_{\perp}-\delta\bm{z}^{-}_{\perp}$. Elsässer variable
increments were determined as
$\delta\bm{z}_{\perp}^{\pm}=\delta\bm{v}_{\perp}\pm\text{sign}(B_{r}^{0})\delta\bm{b}_{\perp}$,
where $B_{r}^{0}$ represents the 30-minute rolling average of the radial
magnetic field component, $B_{r}$, used to determine the polarity of the
background magnetic field. Here, $\bm{z}_{\perp}^{-}$ and $\bm{z}_{\perp}^{+}$
denote inward and outward propagating Alfvén waves, respectively. Magnetic
field data were downsampled, following the application of a low-pass
Butterworth (1930) filter to mitigate aliasing, to match the temporal
resolution of the ion moment data. Magnetic field data were then normalized to
velocity units using a 1-minute moving average applied to the proton density,
$n_{p}$, time series.

Figure 7a,b illustrates two sets of alignment angles: $\theta_{\perp}^{ub(z)}$
in black, as defined by Equation 10, and $\tilde{\theta}_{\perp}^{ub(z)}$ in
blue, calculated according to Equation 13. These sets are referred to as
$\Theta^{z}$ when discussing Elsässer variables and $\Theta^{ub}$ in the
context of velocity-magnetic field fluctuations, with the angle range confined
to $0^{\circ}-90^{\circ}$, consistent with Podesta et al. (2009). The inset of
Figure 7b illustrates the normalized residual energy, $\sigma_{r}$, in red,
and the normalized cross helicity, $\sigma_{c}$, in gray.

At the energy injection range, $\lambda\gtrsim 2\times 10^{4}d_{i}$, a trend
towards tighter alignment at smaller scales is observed, predominantly in
$\Theta^{ub}$. This is accompanied by a monotonic increase in $\sigma_{c}$ and
a shift of $\sigma_{r}$ towards more negative values. These trends are more
pronounced over longer observational intervals, although such extended periods
are beyond the scope of this analysis.

Figure 6: The scale-dependent kurtosis of the magnetic field, denoted as
$K_{5}$ and $K_{2}$, is estimated using 5-point structure functions
($SF_{5}^{n}$) and 2-point structure functions ($SF_{2}^{n}$), respectively.
These are plotted in red, gray and light-blue, respectively, as a function of
scale in units of ion inertial length ($d_{i}$) for the parallel (a),
perpendicular (b), and displacement (c) directions. Power-law fits have been
applied to $K_{5}$ over the region $R_{2}$.

Within the $R_{2}$ range, $\theta^{z}$ exhibits negligible variation with
scale. Conversely, $\tilde{\theta}^{z}$ reveals subtle signatures of enhanced
alignment at $\lambda\lesssim 2\times 10^{3}d_{i}$, coinciding with
$\sigma_{r}$ transitioning from negative to positive values. Simultaneously,
$\Theta^{ub}$ steadily increases, indicating progressive misalignment between
the magnetic field and velocity increments at inertial scales.

The disparity between the two alignment definitions becomes more evident in
the $R_{1}$ range, where $\tilde{\theta}^{z}\propto\lambda^{0.11}$. However,
observations concerning this range should be approached cautiously due to
potential instrumental noise, a matter further explored in Section 6.

What appears to solidly be the case, however, is that an inverse relationship
holds between the alignment angle and the intensity of the field gradients.
This relationship is depicted in the inset panel of Figure 7a, where
$\theta^{ub}$ is plotted across various scale-dependent percentile bins of the
Partial Variance of Increments (PVI) diagnostic,
$\mathcal{I}_{\bm{B}}(t,\ell)=|\delta\bm{B}(t,\ell)|/\sigma_{\bm{B}}$, with
$\sigma_{\bm{B}}$ representing the standard deviation calculated over a moving
window of 1 hour, (Greco et al., 2018). Specifically, at $\lambda\gtrsim
2\times 10^{2}d_{i}$, higher-percentile $P(I_{\bm{B}})$ bins are characterized
by lower average $\theta^{ub}_{\perp}$ values. Similar results were obtained
when considering $\theta^{z}_{\perp}$ and segregating alignment angles based
on the percentile bins of the PVI diagnostic applied to the
$\bm{z^{+}_{\perp}}$ time-series, $I_{\bm{z}^{+}}$.

Under the assumption that the cascade is local in $\lambda$, we investigated
the scale-dependence of the nonlinearity parameter $\chi^{\pm}$, according to
the formulations by B06 and CSM15:
$\chi^{\pm}=(\ell^{\pm}_{||,\lambda}/\lambda^{\pm})/(\delta
z_{\lambda}^{\mp}/V_{a})sin\theta^{z}$. The analysis of outgoing ($\chi^{+}$)
and ingoing ($\chi^{-}$) waves, depicted in gray and red respectively in panel
(c) of Figure 7, reveals that both cascades start weak, with $\chi^{\pm}<1$.
As the cascade progresses towards smaller scales, a significant increase in
the wave to nonlinear times ratio leads to a consistent rise in $\chi^{\pm}$,
continuing until scales nearing the $R_{2}$ range onset.

For the ingoing waves, a transition from weak to strong wave turbulence is
noted, with $\chi^{-}>1$ at $\lambda\approx 3\times 10^{4}d_{i}$. The cascade
remains strong throughout the resolvable portion of the inertial range, with
$\chi^{-}$ being scale-independent, staying close to, yet slightly greater
than, 1.

In contrast, the cascade of outwardly propagating waves remains weak within
the $R_{2}$ range. More specifically, $\chi^{+}$ shows a modest increase from
approximately 0.1 at $\lambda=10^{4}d_{i}$ to around 0.2 at $\lambda=2\times
10^{2}d_{i}$.

The potential inaccuracies in velocity measurements, exacerbated at smaller
scales, along with the limited resolution of velocity field data, caution
against drawing definitive conclusions about the cascades’ nature in the
$R_{1}$ range.

Lastly, it’s noteworthy that the definition of $\chi^{\pm}$, as proposed by
GS95, was also considered. This analysis revealed a scale dependence similar
to that of $\chi^{\pm}$, but with both $\chi^{\pm}$ values being approximately
twice as high.

##  6 Discussion

Recent in-situ observations indicate that the regime canonically identified as
the inertial range comprises two sub-inertial segments, exhibiting distinct
scaling behaviors (Wicks et al., 2011; Chhiber et al., 2021; Sioulas et al.,
2022; Telloni, 2022; Wu et al., 2022; Sioulas et al., 2023; Sorriso-Valvo et
al., 2023).

Building on this insight, our study investigates the anisotropic properties
and higher-order statistics of the two sub-inertial ranges, utilizing a
physically motivated, locally defined coordinate system. Concurrently, we
focus on evaluating the predictions of homogeneous MHD turbulence models,
grounded in the principles of $\it{critical~{}balance}$ and
$\it{dynamic~{}alignment}$, as proposed by CSM15 and MS17.

In the ensuing section, we embark on a detailed comparison with previous
theoretical, observational, and numerical results that contextualizes our
findings.

###  6.1 Investigating the Impact of Imbalance and Expansion on the Higher-
Order Statistics

Phenomenological treatment of homogeneous MHD turbulence (e.g., GS95, B06,
CSM15, MS17) is usually performed under the simplifying assumption of
negligible cross helicity. However, the statistical properties of solar wind
turbulence vary significantly with the degree of Elsässer and Alfvénic
imbalance (Podesta & Borovsky, 2010; Chen et al., 2013; Wicks et al., 2011,
2013b; Bowen et al., 2018; Andrés et al., 2019; Sorriso-Valvo et al., 2021;
Sioulas et al., 2023; D’Amicis et al., 2022; McIntyre et al., 2023). Various
models have been proposed as modifications to the frameworks of GS95 and B06,
incorporating different assumptions about the turbulent cascade to address the
imbalance in oppositely directed Alfvénic fluxes (Lithwick et al., 2007;
Beresnyak & Lazarian, 2008; Chandran, 2008; Perez & Boldyrev, 2009; Podesta &
Bhattacharjee, 2010; Schekochihin, 2022). As these models omit considerations
of intermittency, they will not be elaborated upon in the ensuing discussion.

Figure 7: Alignment angles (a) $\Theta^{ub}(\lambda)$ and (b)
$\Theta^{z}(\lambda)$, constrained to the range of
$0^{\circ}\leq\Theta^{ub(z)}\leq 90^{\circ}$. The black curves depict
alignment angles calculated with Equation 10, whereas the blue curves are
based on Equation 13. The inset in panel (a) highlights the inverse
relationship between alignment angle and field gradients by color-coding
$\theta^{ub}$ across scale-dependent percentile bins of $I_{\bm{B}}$. The
$i$-th bin at scale $\lambda$ is defined as $[10(i-1),10i)$, for $i=1,...,10$.
The inset in panel (b) displays $\sigma_{r}(\lambda)$, in red, and
$\sigma_{c}(\lambda)$, in gray. Panel (c) shows the nonlinearity parameter,
$\chi^{\pm}(\lambda)$, for outwardly ($\delta z^{+}$) and inwardly ($\delta
z^{-}$) propagating waves, depicted with gray and red lines respectively. For
all cases, results were derived by first calculating values for each interval
independently and then computing a scale-dependent median.

From an observational standpoint, the extent to which imbalance impacts
higher-order statistics in MHD turbulence, is largely unexplored territory.
Previous investigations have predominantly focused on categorizing findings
based on wind speed or analyzing mixed Fast/Slow streams (Horbury & Balogh,
1997b; Mangeney, 2001; Salem et al., 2009; Chhiber et al., 2021; Wu et al.,
2023), with a recent shift in interest towards magnetic compressibility
(Palacios et al., 2022). However, these studies often do not explicitly detail
the degree of Elsässer imbalance in the dataset, making direct comparisons
with our observations challenging. Nonetheless, it is worth pointing out that
the scaling profiles observed in the $R_{2}$ range, especially for the
parallel and perpendicular components, qualitatively align with results
presented in Osman et al. (2012).

In terms of numerical simulations, both homogeneous and inhomogeneous setups
have been employed to examine the influence of imbalance on MHD turbulence
statistics. Studies have primarily concentrated on spectral properties,
revealing inertial range scalings of $\alpha_{\lambda}=-3/2$ (Perez et al.,
2012; Chandran & Perez, 2019; Grappin et al., 2022; Meyrand et al., 2023).

Shi et al. (2023) investigated the effects of Elsässer imbalance on higher-
order statistics using both homogeneous and EBM simulations, each initialized
with varying degrees of imbalance. It was shown that EBM simulations
initialized with $|\sigma_{c}|\approx 1$ result in higher-order moment scaling
exponents consistent with the predictions of the CSM15 model. Conversely,
simulations with lower values of $\sigma_{c}$ exhibited scaling exponents
showing a linear dependence on order $n$. The study also highlighted
significant differences in higher-order scaling exponents between homogeneous
and EBM simulations, even when initialized with identical levels of imbalance.
For example, scaling exponents in homogeneous runs with $\sigma_{c}\approx 0$
were found to be convex functions of order and closely resembled the CSM15
model, in line with previous investigations (Chandran et al., 2015; Mallet et
al., 2015; Palacios et al., 2022). Additionally, Shi et al. (2023) noted that
in EBM simulations, the scaling properties displayed variations when higher-
order moments were calculated from increments sampled in directions other than
radial.

The latter observation is consistent with the findings of Verdini & Grappin
(2015), who conducted a comparative analysis of balanced homogeneous and EBM
simulations. They observed that while the homogeneous simulations displayed
three-dimensional anisotropy, in agreement with B06—a finding further
corroborated by Mallet et al. (2016)—the EBM simulations demonstrated
axisymmetry relative to $\bm{B}_{\ell}$ and did not exhibit three distinct
inertial range scaling laws. Specifically, EBM simulations with increments
measured along the radial direction demonstrated spectral scalings of
$\alpha_{\lambda}\approx\alpha_{\xi}\approx-3/2$ in both the perpendicular and
displacement components. However, measurements in the transverse direction
revealed scalings of $\alpha_{\lambda}\approx\alpha_{\xi}\approx-5/3$. In both
cases, the parallel component lacked convincing scaling properties, although
it exhibited a slightly steeper spectrum compared to the perpendicular
components.

These findings lead to the following interpretation: The large-scale flow,
being radial rather than uniform, cannot be negated by a Galilean
transformation. Consequently, the expansion introduces an additional axis of
symmetry and fosters a scale-dependent competition between the mean-field and
radial axes (Völk & Aplers, 1973). Intuitively the effects of the expansion
should be important when the non-linear time, $\tau_{nl}$ is slower than the
expansion of the solar wind, $\tau_{exp}=R/V_{SW}$, $\tau_{exp}\leq\tau_{nl}$.
Given that $\tau_{nl}$ is scale-dependent, it logically follows that the
effects of expansion are more pronounced at larger scales777Note, however,
that observational evidence suggests that expansion can contaminate the
turbulence statistics even within inertial scales (Verdini et al., 2018,
2019). From this discussion, it becomes evident that the expansion has the
potential to modify the local 3D anisotropy in a scale-dependent manner.

At the resolution currently achievable in (R)MHD simulations, meaningful
comparisons are possible with the larger-scale end of the inertial range,
$R_{2}$, spanning $200-6000d_{i}$. Focusing on $\zeta_{\lambda}$ and setting
aside anisotropy, our findings are consistent with Shi et al. (2023)’s
imbalanced EBM results and, consequently, with the models of CSM15 and MS17.
Nevertheless, the observed discrepancies in the scaling exponents of the
parallel and displacement components within $R_{2}$ could hint at the
influence of expansion effects, suggesting a scenario where the dominant axis
of symmetry is a mix of both $B_{\ell}$ and the radial axis, with the
contribution of each being scale-dependent.

Shifting focus to the scaling exponents in $R_{1}$ (spanning $8-100d_{i}$),
recent theoretical work suggests that under conditions of strong imbalance,
generalized helicity conservation may hinder turbulent energy transfer to
kinetic scales (Passot et al., 2018; Passot & Sulem, 2019; Meyrand et al.,
2021; Passot et al., 2022). The “helicity barrier” effect could influence both
spectral and potentially higher-order moment scalings at the smaller end of
the inertial range. Homogeneous hybrid-kinetic simulations, initialized under
strong imbalance conditions to capture this effect (Squire et al., 2022a),
exhibit spectral exponents for the parallel and perpendicular magnetic field
components, $a_{\ell_{||}}\approx-2$ and $a_{\lambda}\approx-5/3$
respectively, in line with observations in the $R_{1}$ range.

While our current data does not definitively link our observations to the
helicity barrier effect, it underscores the necessity for more comprehensive
numerical studies focusing on the higher-order statistics of strongly
imbalanced turbulence.

In conclusion, this discussion underscores the significant impact of imbalance
and expansion effects on MHD turbulence statistics, indicating that the solar
wind might not provide an ideal laboratory for evaluating homogeneous MHD
turbulence models. This realization calls for a prudent application of
homogeneous, balanced turbulence models in the analysis of solar wind
observations and highlights the imperative for more sophisticated theoretical
modeling and refined data interpretation techniques.

###  6.2  Critical Balance (CB)

Using balanced RMHD simulations, Mallet et al. (2015) demonstrated that
although the distributions of $\tau_{nl}$ and $\tau_{A}$ are not self-similar,
their ratio, $\chi$, maintains a scale-invariant distribution within the
inertial range. Chhiber et al. (2020) employed balanced incompressible MHD
simulations to show that while the $\chi$ distribution peaks at $\chi\approx
1$, it is asymmetric and skewed towards $\chi\geq 1$. Further reinforcing the
Chhiber et al. (2020) findings, Oughton & Matthaeus (2020) highlighted that,
despite RMHD simulations producing results claimed to support CB, the
similarity between RMHD’s $\chi\geq 1$ requirement and the CB condition of
$\chi\approx 1$ has led to some confusion in differentiating these two
theoretical frameworks.

Chen (2016) utilized an extensive dataset of fast wind streams with moderate
cross-helicity ($\sigma_{c}\approx 0.6$) from the outer heliosphere to
investigate the scale dependence of the non-linearity parameter. They found
$\chi$ to be scale-independent across the inertial range, maintaining a value
around $\chi\sim 1$. Due to the lower resolution in velocity data and a less
pronounced imbalance compared to our dataset, they assumed identical
statistical properties for the two Elssässer fields, enabling them to estimate
$\chi=(\ell_{||}/\lambda)(\delta b/V_{A})$ using solely magnetic field data.

However, numerical simulations by Beresnyak & Lazarian (2008, 2009) suggest
that with increasing imbalance, statistical properties (i.e., amplitudes,
coherence lengths) of the two Elsässer species diverge progressively.
Considering the strong imbalance in our dataset, we employed a more refined
approach, computing
$\chi^{\pm}=(\ell^{\pm}_{||,\lambda}/\lambda^{\pm})/(\delta
z_{\lambda}^{\mp}/V_{a})sin\theta^{z}$.

Our analysis reveals a strong cascade for the inwardly propagating waves, with
$\chi^{-}$ remaining scale-independent across the inertial range, maintaining
a value slightly above unity. In contrast, the outwardly propagating waves
exhibit a weaker cascade, with $\chi^{+}$ increasing from $0.1$ at
$\lambda\approx 10^{4}d_{i}$ to $0.2$ at $\lambda\approx 10^{2}d_{i}$

At this point, it’s important to recognize two key factors that might affect
the accuracy of our $\chi^{\pm}$ estimates. First, there’s a prevailing
assumption that $z^{\pm}$ structures are primarily sheared by counter-
propagating $z^{\mp}$ wavepackets of similar perpendicular scale, which
implies a cascade that is local in $\lambda$. However, this notion is
challenged by the work of Schekochihin (2022) who put forward a model of
imbalanced turbulence that consists of two strong, “semi-local” cascades: one
local in $\lambda$ for the stronger field and another, local in $\ell_{||},$
for the weaker field. The implications of this model cast doubt on our
estimates of $\chi^{\pm}$ that is predicated on the concept of scale locality.

Furthermore, the transition from weak to strong turbulence, characterized by a
change from $\chi\ll 1$ and $\alpha_{\lambda}\approx-2$ to $\chi\approx 1$ and
$\alpha_{\lambda}\approx-3/2$, is a pivotal aspect of the CB discussion. This
transition, observed in balanced RMHD shell-model simulations Verdini &
Grappin (2012), 3D incompressible MHD simulations (Meyrand et al., 2016), and
recently in the Earth’s magnetosheath (Zhao et al., 2023), remains unreported
in the solar wind 888Even though such a transition was speculated in recent
works (Telloni, 2022; Wu et al., 2022; Sioulas et al., 2023).

Our results indicate a transition in the cascade of the ingoing wave from weak
to strong turbulence, with $\chi^{-}>1$ at $\lambda\approx 3\times
10^{4}d_{i}$, yet without capturing the anticipated scaling transition. We
speculate, below, that this could be related to the effects of “anomalous
coherence”. More specifically, in the context of homogeneous MHD, nonlinear
interactions between counterpropagating waves are uncorrelated and transient,
limited to the duration of encounters. In the solar wind, however, nonlinear
dynamics are complicated by anomalous coherence, a phenomenon arising from
non-WKB reflection of outwardly propagating fluctuations (Velli et al., 1989,
1990; Hollweg & Isenberg, 2007). A key aspect of this effect is the presence
of an “anomalous” reflected component, $z^{-}_{a}$, in addition to the
“classical” component $z^{-}_{c}$, which remains stationary relative to the
$z^{+}$ frame, coherently shearing it throughout its lifetime (Verdini et al.,
2009; Chandran & Perez, 2019). Under strong imbalance and inhomogeneity,
$z^{-}_{a}$ can assume a leading role at large scales, altering the
phenomenology of the energy cascade and leading to distinct spectral
behaviors: a $1/f$ scaling for outwardly propagating waves and an $f^{-3/2}$
scaling for inwardly propagating waves (Velli et al., 1989; Perez & Chandran,
2013; Meyrand et al., 2023).

In summary, our results suggest that at inertial scales, outgoing waves
experience a weak cascade, while ingoing waves undergo a strong one, closely
resembling the CB condition ($\chi^{-}_{\xi}\approx 1$). However, given the
complexities previously discussed and the uncertainties inherent in our
measurements, we advise interpreting these findings with a degree of caution.

###  6.3  Scale Dependent Dynamic Alignment (SDDA)

Several numerical investigations into homogeneous (R)MHD (Mason et al., 2006;
Perez et al., 2012, 2014; Cerri et al., 2022) have provided ample evidence for
alignment signatures spanning a significant portion of the inertial range.
However, Beresnyak (2012, 2014) have suggested that this observed alignment
increase may be a finite-range phenomenon closely linked to dynamics at the
outer scale.

Observational studies using data sampled at 1 AU have provided evidence of
alignment at large, energy-containing scales. However, it has been observed
that this trend towards increasing alignment diminishes at inertial scales
(Podesta, 2009; Hnat et al., 2011; Wicks et al., 2013a; Parashar et al.,
2019). This trend persists even in data intervals specifically chosen to
mitigate the effects of solar wind expansion (Verdini et al., 2018). It has
been noted, however, that small errors in velocity vector measurements, due to
instrumental limitations, can lead to significant errors in alignment angle
measurements, even at large scales (Podesta et al., 2009).

Recently, Parashar et al. (2020) explored the scale dependence of several
alfvénicity diagnostics during E1 of PSP. Their findings suggest that
$\sigma_{c}$ starts decreasing, with $\sigma_{r}$ increasing, at scales
considerably larger than those observed at 1 AU (Podesta & Borovsky, 2010),
despite high alfvénicity at large scales. This observation aligns with HELIOS
observations (Tu et al., 1990) and has been attributed to the substantial
energy found in velocity shears in the inner heliosphere (Ruffolo et al.,
2020). Specifically, shear disrupts an initial spectrum of high cross helicity
by injecting equal amounts of the two Elsässer energies (Roberts et al., 1987;
Goldstein et al., 1989; Roberts et al., 1992).

The scale-dependence of the alignment angles correlates directly with that of
$\sigma_{c}$ and $\sigma_{r}$, as only two out of these four quantities are
independent. Specifically, the formal relationship between imbalance, residual
energy, and alignment, as described by
$cos\theta^{z}_{\perp}=\sigma_{r}/(1-\sigma_{c}^{2})^{1/2}$ and
$cos\theta^{ub}_{\perp}=\sigma_{c}/(1-\sigma_{r}^{2})^{1/2}$ , indicates that
the development of both Elsässer imbalance and residual energy, i.e., a
monotonic increase in $|\sigma_{c}|$ and $|\sigma_{r}|$ towards smaller
scales, is necessary for SDDA to emerge towards smaller scales (Wicks et al.,
2013a; Schekochihin, 2022).

Such trend is evident in Figure 7, where alignment signatures become apparent
only when $\sigma_{c}$ exhibits a monotonic increase at scales
$\lambda<10^{4}d_{i}$, or when $\sigma_{r}$ becomes positive, leading to a
monotonic increase in $|\sigma_{r}|$ at scales $\lambda<2\times 10^{3}d_{i}$.
Although instrumental noise might influence the latter trend, as discussed in
Bourouaine et al. (2020), the observed behavior at energy injection scales
aligns with 1 AU observations (Wicks et al., 2013a). In contrast, in the
inertial range, $\theta^{z}(\approx 35^{\circ})$ remains roughly scale-
independent.

These observations raise a critical question: Is the observed scale-dependence
of the alignment a reflection of actual physical processes, or might it simply
be a consequence of instrument characteristics?

Potential physical mechanisms may encompass interactions between compressive
and non-compressive modes (Cho & Lazarian, 2003; Chandran, 2005), ideal MHD
instabilities manifesting in the solar wind, including the Kelvin-Helmholtz
instability (Malagoli et al., 1996), the cessation of the aligning cascade due
to the tearing instability (Mallet & Schekochihin, 2017; Boldyrev & Loureiro,
2017; Comisso et al., 2018), or even the solar wind’s inherent inhomogeneity
resulting in non-WKB reflections and a reduction in cross-helicity.

The extent to which instrumental noise influences these observations remains a
crucial, yet unresolved, concern, emphasizing the necessity for a careful
interpretation of observational data. While our analysis cannot definitively
assert the nature of SDDA at small scales, it provides compelling
observational evidence suggesting an inverse relationship between alignment
angles and the intensity of field gradients, thereby corroborating the
numerical results of Mallet et al. (2015). As such, our observations lend
support to the CSM15 model, suggesting that the physical basis of alignment
lies in the mutual shearing of Elsässer fields during imbalanced collisions
between counterpropagating wave packets.

Before concluding, another topic related to our findings deserves further
discussion. Our observations indicate strong anti-alignment between Elsässer
species at larger scales, with $\sigma_{c}\approx 1$ and typically small and
negative $\sigma_{r}$, implying, as suggested by Wicks et al. (2013a),
$cos\theta^{z}_{\perp}\approx 180^{\circ}$ at $\lambda\approx 10^{4}d_{i}$. It
is natural, then, to inquire how this picture would be modified in the case of
globally balanced streams. For instance, as emphasized in CSM15, outer-scale
fluctuations in the context of globally balanced turbulence are not expected
to be strongly aligned. This suggests that while imbalanced turbulence may
exhibit alignment saturation at larger scales, balanced turbulence has the
potential for further alignment, assuming that $\sigma_{c}$ has room to
increase at smaller scales. Thus, the dynamics of alignment across scales in
balanced turbulence present an intriguing area for future research,
particularly in comparing the extent of this alignment process with that in
imbalanced streams. This aspect warrants further investigation and will be the
focus of future work.

###  6.4 Exploring the Efficacy of the CSM15 & MS17 Models in the Context of
Imbalanced Turbulence

Our results indicate that the intermittent scalings of the $n$th-order
conditional structure functions in the direction perpendicular to both the
local mean field and the fluctuation directions closely align with the
theoretical frameworks proposed by CSM15 and MS17. However, our dataset
exhibits notable characteristics that diverge from the assumptions underlying
these models. These include strong imbalance in the fluxes of the Elsässer
species and a prevalence of eddies conforming to a field-aligned tube
topology, which, on average, do not display increasing alignment towards
smaller scales. At first glance, these observations might seem contradictory
to the expectations set forth in the aforementioned models. Nevertheless, upon
closer examination of the model’s fundamental assumptions, it becomes apparent
that these empirical findings are not inconsistent with the model’s
theoretical framework.

In Section 2, we examine the foundational assumptions of the CSM15 and MS17
models. These assumptions encompass (1) negligible cross helicity at energy
injection scales and (2) the formation of eddies following a current sheet
topology, characterized by a volume filling factor $f_{cs}\propto\lambda$,
alongside the dynamic alignment of vector field fluctuations at smaller
scales. The models incorporate alignment as an intermittency effect, resulting
in (1) the inherent introduction of local Elsässer imbalance and (2) the lack
of a requirement for the “average” eddy to adhere to a 3D anisotropic current
sheet topology. The models propose that it is the eddies in the tails of the
PDFs that are expected to exhibit this topology, showing increased alignment
at smaller scales. In essence, an inverse correlation between alignment and
the intensity of field gradients is anticipated at any given scale.

When second-order moments are employed to examine the statistical shape of
eddies, the high-amplitude, current-sheet-like structures found in the tails
of the increment PDFs are typically obscured during the averaging process by
more common, lower-amplitude, tube-like eddies. However, as higher moment
orders are considered, these high-amplitude structures gain more prominence in
the ensemble average, thereby significantly impacting the profile of the
scaling exponents.

Given that fluctuations in the perpendicular component are energetically
dominant, any “side effects” stemming from the disruption of anisotropy due to
expansion, and possibly from constraints imposed by spherical polarization
(Matteini et al., 2024), do not seem to significantly impact the scaling
characteristics of the perpendicular component. However, they do appear to
influence the scaling profiles of the displacement and parallel components.
Hence, incorporating these effects into the existing homogeneous models could
potentially enhance the scaling predictions of these models and offer valuable
insights into the nature of intermittent turbulence in the expanding solar
wind.

###  6.5 Observations of the outer scale

At large energy-containing scales, as depicted in Figure 3, the parallel and
perpendicular components display rough equipartition in fluctuating energy.
However, the displacement component appears somewhat more energetic in
comparison. Consequently, eddies deviate from isotropy, exhibiting a subtle
compression along the fluctuation direction.

Our findings contrast with those of Chen et al. (2012), who observed eddy
elongation along the displacement direction ($\hat{\xi}$) at large scales in
fast solar wind data at 1.4 AU.

To delve deeper into these findings, Verdini & Grappin (2015) conducted a
comparative analysis using 3D MHD homogeneous and EBM simulations. They found
isotropic eddies at energy injection scales in non-expanding simulations.
However, in expanding simulations where increments were measured along the
radial direction (typical for single-spacecraft missions at 1 AU), results
aligned with those reported by Chen et al. (2012). Interestingly, large-scale
anisotropy disappeared, and eddies appeared isotropic when increments were
measured in non-radial directions. This was interpreted as an effect of
expansion, which preferentially dampens the radial component of magnetic field
fluctuations relative to the azimuthal ones, confining fluctuations primarily
to the plane orthogonal to the radial direction and leading to an anisotropic
energy distribution among the field components. This phenomenon was observed
in 3D EBM simulations (see e.g., Dong et al., 2014).

In a more recent study, Verdini et al. (2018) analyzed a decade of data from
the Wind spacecraft and identified a correlation between large-scale
anisotropy and variance anisotropy, defined as
$E=(b^{2}_{T}+b^{2}_{N})/b^{2}_{R}$, where $b$ represents the rms amplitude of
fluctuations. They found that intervals corresponding to the “strong”
expansion dataset (defined by $E>2$) exhibited eddy elongation along the
displacement direction, consistent with Chen et al. (2012). Conversely,
intervals from the “weak” expansion dataset (with $E\leq 2$) showed eddy
elongation along the perpendicular direction ($\hat{\lambda}$).

Based on these results, Verdini et al. (2018) suggested that PSP, due to its
unique orbit– allowing measurements perpendicular to the radial direction
during its near-sun phase– would detect isotropic eddies at energy injection
scales.

To overcome the limitations of single spacecraft measurements, Vech & Chen
(2016) adopted a multispacecraft approach, enabling the separation of
measurements along both radial and transverse directions. This strategy
facilitated the isolation of expansion, solenoidality, and the mean magnetic
field effects. Their investigation underscored the dominant role of the
solenoidality constraint (Turner et al., 2011) over expansion in contributing
to the observed variance anisotropy, E (see e.g., Horbury & Balogh, 2001).
They further noted that while some anisotropy, as observed with radial
increments, stemmed from expansion, a reduced yet noticeable degree of
anisotropy persisted when measurements were conducted along the transverse
direction. This suggests the involvement of additional mechanisms in shaping
the large-scale 3D anisotropy.

Recent in-situ observations have shed light on the decay observed in the
radial component within the inner heliosphere, indicating that it cannot be
solely attributed to expansion effects. Specifically, analysis of data from
PSP and HELIOS by Tenerani et al. (2021) demonstrated that the rms of
fluctuations in the radial component decays at a slower rate compared to that
of the perpendicular component. This phenomenon was further clarified by
Matteini et al. (2024), who attributed it to the tendency of magnetic field
fluctuations in the solar wind to evolve towards a state of spherical
polarization. The spherical polarization imposes constraints on the radial
component’s rms fluctuations, leading to a decay described, particularly at
large scales, by $b_{r}\sim\bm{b}/2B$ (see also, Squire et al., 2020; Mallet
et al., 2021).

Returning to our findings, we observe only a slight deviation from isotropy at
the outer scale, which is notable considering that our dataset falls within
the strong expansion category, with an average variance anisotropy of
$E\approx 4.26\pm 3.51$. Additionally, the sampling angle consistently falls
within the range of $\theta_{R}\in[160^{\circ},180^{\circ}]$, indicating
quasi-radial sampling. Despite identifying several intervals among multiple
PSP encounters that exhibit isotropic large-scale eddies, no clear correlation
has emerged from this preliminary analysis between prevailing plasma
signatures (e.g., $E$, $\sigma_{c}$, $\sigma_{r}$, $\theta_{R}$) and the
occurrence of such intervals.

Taking into account the subtleties revealed by the recent observations
discussed above, we must recognize our current inability to offer a
satisfactory explanation for the observed configuration of the large-scale
eddies in our dataset. However, the differences noted in the near-Earth and
near-Sun eddies could indicate preferential damping in the fluctuations of the
displacement component of the magnetic field. To clarify this aspect, it would
be worthwhile to explore the development of the large-scale eddies at varying
heliocentric distances.

###  6.6 Isotropization of eddies at small scales

Figure 4 illustrates that within the $R_{1}$ range, the eddies display
increasing anisotropy, resembling ribbon-like structures towards smaller
scales. However, the trend of increasing aspect ratio ceases at
$\lambda\approx 2d_{i}$, at which point the eddies transition toward a quasi-
isotropic state. In the following, we discuss two potential mechanisms that
could explain these observations.

The observed transition towards isotropy at smaller scales is consistent with
the idea that thin, long-lived current sheets generated by the turbulent
cascade can be disrupted by the tearing instability and subsequent
reconnection (Furth et al., 1963). Specifically, when the maximum growth rate
of the Coppi et al. (1976) mode, $\gamma_{t}$, representing the fastest
tearing mode in an MHD sheet, becomes comparable to the non-linear cascade
time $\tau_{nl}$ ($\gamma_{t}\tau_{nl}\gtrsim 1$), the stability of the
current sheets is compromised (Pucci & Velli, 2014; Uzdensky & Loureiro,
2016). The instability leads to the fragmentation of the dynamically forming
sheets into flux ropes, which exhibit isotropy in the plane perpendicular to
the magnetic field. This phenomenon is expected to occur at scale,
$\lambda_{D}$, beyond which the nature of the MHD cascade undergoes a
significant transformation. The disruption of the current sheets affects the
dynamic alignment, accelerating the turbulent cascade and resulting in a
noticeable steepening of the power spectrum (Mallet et al., 2017; Loureiro &
Boldyrev, 2017; Boldyrev & Loureiro, 2017).

While this mechanism appears feasible for balanced turbulence, as evidenced by
observations in the solar wind (Vech et al., 2018) and more recent findings in
3D fully-compressible (Dong et al., 2022) and reduced (Cerri et al., 2022) MHD
simulations, it’s important to acknowledge that the cutoff of the inertial
range in imbalanced turbulence might also be influenced by kinetic effects.
Therefore, despite the observation of a sub-ion-scale range mediated by
magnetic reconnection in 2.5D hybrid-kinetic simulations (Cerri & Califano,
2017; Franci et al., 2017), other kinetic-scale mechanisms in imbalanced
turbulence could potentially contribute to eddy isotropization. For instance,
an alternative explanation for the isotropization of the eddies at small
scales could be associated with the recently discovered “Helicity Barrier”
mechanism (Passot et al., 2018; Passot & Sulem, 2019; Meyrand et al., 2021;
Passot et al., 2022).

As discussed in the introduction, strongly magnetized (low-beta) collisionless
plasmas exhibit nonlinear conservation of both energy and cross helicity.
However, the conserved quantity in reality is termed “generalized helicity.”
At $k_{\perp}\rho_{i}\lesssim 1$, this corresponds to the cross helicity
following a forward cascade, which conservatively transforms into magnetic
helicity at $k_{\perp}\rho_{i}\gtrsim 1$, undergoing an inverse cascade.
Consequently, an imbalanced cascade arriving from the inertial range faces a
complication—the sudden need to reverse the direction of the generalized
helicity cascade. The helicity barrier impedes the viability of a constant-
flux cascade, leading to an accumulation of energy in the stronger Elsässer
field. This accumulation shortens $\tau_{nl}$, reducing the parallel
correlation length, in line with the CB theory, to the extent that turbulent
energy is redirected into an (ICW; (see, Stix, 1992)) spectrum. This mechanism
opens up a new dissipation channel via the ion-cyclotron resonance (Squire et
al., 2022b).

ICWs are commonly observed in the nascent solar wind, particularly during
intervals marked by (anti)alignment between the mean magnetic field and solar
wind flow direction (Bowen et al., 2020b). A strong correlation exists between
the presence of ion-scale waves and the level of imbalance of fluctuations at
inertial scales (Zhao et al., 2022; Bowen et al., 2023). The presence of ICWs
can significantly impact the power spectra of magnetic fields at ion kinetic
scales (Bowen et al., 2020b; Shankarappa et al., 2023). Specifically, the bump
observed in the parallel spectrum just before the transition region has been
attributed to the presence of ICWs, suggesting that the isotropization of the
eddies could be a consequence of the helicity barrier mechanism.

The main emphasis of this analysis is on the inertial and energy injection
scales. The extended interval size allows for more reliable estimates of
second-order moments and, consequently, the anisotropic curves presented in
Figure 3. However, no effort has been made to account for the energy
contribution ICWs. Due to the strong correlation between the presence of ion-
scale waves and Elsässer imbalance, simply disregarding intervals with ICW
wave signatures (Duan et al., 2021; Zhang et al., 2022) may impede our
investigation into kinetic-scale turbulence statistics in strongly imbalanced
intervals. An alternative approach would entail identifying and eliminating
the energy attributed to ICWs from the observed energy spectrum (see e.g.,
Shankarappa et al., 2023; Wang et al., 2023), and examining the resulting
anisotropy based on the parallel and perpendicular spectra. This will be the
focus of an upcoming study.

###  6.7 Can the trace PSD be interpreted as the perpendicular PSD ?

The theoretical models discussed in Section 1 provide scaling predictions for
the parallel, perpendicular, and displacement components of fluctuating
fields.

However, the angle between the solar wind flow and the magnetic field, as
observed by spacecraft, can significantly influence whether fluctuations in
measured quantities vary parallel or perpendicular to the magnetic field. Due
to the Parker Spiral configuration (Parker & Tidman, 1958), the baseline value
of the angle between the solar wind flow and the magnetic field, denoted as
$\Theta_{VB}$, increases with heliocentric distance. Consequently, spacecraft
measuring magnetic field fluctuations at 1 AU are more likely to detect
fluctuations perpendicular to the mean magnetic field direction.

Considering the strongly anisotropic nature of the turbulent cascade, with the
majority of power associated with perpendicular wavenumbers (Shebalin et al.,
1983; Montgomery & Turner, 1981; Horbury et al., 2008), observational studies
have traditionally estimated the trace PSD and interpreted these measurements
as representative of the perpendicular spectrum. However, as the PSP moves
closer to the Sun, both the flow and the magnetic field become predominantly
radial. Consequently, PSP often detects variations parallel -resulting in a
deficit of measurements perpendicular- to the magnetic field. This can impact
the statistical signatures of MHD turbulence, including intermittency (Sioulas
et al., 2022), estimates of correlation lengths (Cuesta et al., 2022), etc. As
shown in Figure 2, the deficit caused by the sample size of the perpendicular
fluctuations due to quasi-parallel sampling can also result in strong
deviations between the perpendicular and trace PSD. This effect becomes even
more important at later PSP encounters, and thus caution should be exercised
when trying to utilize the trace PSD to compare with perpendicular PSD
predictions of theoretical models.

###  6.8 Assessing 2-Point Structure Functions for Small-Scale Turbulence
Analysis

The departure from monofractal statistics, exemplified by the adoption of
$SF^{n}_{5}$, underscores the inadequacy of the $SF^{n}_{2}$ method for
statistical analysis of MHD turbulence at smaller spatial scales.
Specifically, the $SF^{n}_{2}$ method lacks accuracy in capturing the scaling
behavior at these scales, where steep scaling in the power spectra and higher-
order moments are frequently seen. This deficiency can result in imprecise
estimations of all associated intermittency metrics, leading to potentially
erroneous interpretations of the nature of the turbulent cascade (see also 3D
kinetic simulations by Cerri et al., 2019).

The prolonged duration of the intervals examined in this study, which may not
be optimal for a concentrated analysis at kinetic scales, in conjunction with
the limitations posed by the presence of ICWs as discussed in Section 5.2,
underscores the need for further investigation. Future studies should consider
employing either 5-point or wavelet-derived structure functions, with a
specific emphasis on distinguishing between balanced and imbalanced turbulence
streams (see, e.g., Bowen et al., 2023). This comprehensive exploration is
essential for achieving a deeper and more precise comprehension of the fractal
properties of MHD turbulence at kinetic scales.

###  6.9 Limitations

In addition to the limitations associated with velocity measurements discussed
earlier, it is essential to acknowledge further inherent limitations in the
analysis presented in this study.

####  6.9.1 Finite sample size effects

The proper study of MHD turbulence hinges on the ability to sample plasma from
a common solar source, typically a single solar wind stream, and gather a
sufficiently large sample size for statistical analysis. Spectral properties
alone are insufficient for assessing scale invariance and fractal properties;
higher-order moments are necessary. While evaluating structure functions is
generally straightforward, estimating scaling exponents presents pitfalls. The
primary concern arises from increased sensitivity to rare and large events as
the order, p, increases. This can lead to finite sample effects dominating the
analysis, especially as emphasis shifts to the poorly sampled tails of the
distribution with higher orders. Consequently, higher-order moments become
susceptible to outliers, rendering estimates of scaling exponents increasingly
unreliable (Dudok de Wit et al., 2013; Palacios et al., 2022). As a rule of
thumb, it is generally deemed safe to compute structure functions up to a
certain order, typically defined as $p_{\text{max}}=\log N-1$, where N
represents the sample size (Dudok de Wit et al., 2013).

In our analysis, these challenges are further compounded by two factors.
Firstly, we employ conditional analysis to estimate higher-order moments in
three physically motivated directions, resulting in the exclusion of a
significant portion of increments falling outside specified angle ranges.
Secondly, we utilize 5-point structure functions, where the way increments are
taken leads to a larger portion of the time series being discarded due to edge
effects.

Moreover, PSP data introduces added complexity compared to 1AU measurements,
as the distance from the Sun rapidly changes, causing large variations in the
rms of fluctuations between intervals sampled at different heliocentric
distances. Consequently, very long intervals cannot be utilized effectively.

Despite these limitations, the present analysis could be significantly
improved by adopting the method described in Palacios et al. (2022), where a
large sample of increments from non-contiguous solar wind streams with similar
characteristics can be utilized to construct the PDFs needed to obtain higher-
order moments. This approach will be the focus of future work.

####  6.9.2 Scaling exponents of Elsässer fields.

In this study, our primary focus has been on the higher-order moments derived
from the magnetic field timeseries. However, it’s crucial to recognize that
the fundamental variables in MHD are the Elsässer fields—rather than $\bm{B}$
and $\bm{V}$—due to their conservatively cascading energies. Indeed, the
scaling predictions provided by the CSM15 and MS17 models pertain to the
scaling exponents of the Elsässer field increment moments. Therefore, a more
direct comparison with these models would entail estimating the moments of
increments in $\bm{z}^{\pm}$ (see e.g., Palacios et al., 2022). However,
adopting this approach would require downsampling the magnetic field
timeseries to synchronize with the cadence of the velocity field data. This
would lead to a notable reduction in sample size, and render the estimation of
anisotropic higher-order moments unfeasible with our current dataset.

Nevertheless, as discussed in CSM15, the regions contributing dominantly to
both types of structure functions are those where $\delta\bm{z}^{\pm}$
exhibits exceptional magnitudes. In these regions, one Elsässer fluctuation,
e.g., $\delta\bm{z}^{+}$, typically dominates over the other, leading to
$\delta\bm{b}\approx(1/2)\delta\bm{z}^{+}$. Therefore, given the significant
imbalance in our dataset, the scaling exponents estimated for the magnetic
field timeseries can provide a reasonable approximation for the scaling
exponents of the dominant (outgoing) Elsässer field.

####  6.9.3 Switchbacks

The near-Sun solar wind environment is characterized by the prevalent
occurrence of Switchbacks, a subset of predominantly Alfvénic fluctuations
with amplitudes significant enough to cause the magnetic field to reverse its
direction abruptly, resulting in a local field polarity reversal and a
corresponding radial velocity jet (Matteini et al., 2014; Horbury et al.,
2018; Bale et al., 2019).

The question arises as to what extent these sudden reversals impact our
ability to accurately estimate the local magnetic field, and consequently, the
scaling exponents of the parallel and displacement components.

While excluding switchbacks from the analysis could potentially address this
concern, it’s noteworthy that the majority of our samples for the two
perpendicular components originate from substantial kinks in the magnetic
field time series, as these events lead to large $\Theta_{BV}$ angles.
Therefore, no attempt has been made to further clarify this aspect. However,
it is reassuring to note that the scaling exponent profiles obtained for a
substantial dataset of imbalanced Wind observations at 1 AU, where significant
kinks in the magnetic field time series typically diminish and
switchbacks/switchback patches transition into microstreams (Horbury et al.,
2023; Soni et al., 2024), are qualitatively consistent with those reported in
the current analysis.

##  7 Conclusions and Summary

We analyzed in-situ observations from a highly Alfvénic stream captured during
Parker Solar Probe’s first perihelion to assess the predictions of MHD
turbulence models grounded on the principles of “Critical Balance” and “Scale-
Dependent Dynamic Alignment”. Our objective was to assess the extent to which
the conjectures made and predictions derived by these models align with in-
situ solar wind observations and establish solid observational benchmarks for
the testing and refinement of MHD turbulence phenomenologies.

The main findings of our study can be summarized as follows:

At the outer scale, $\lambda\gtrsim 2\times 10^{4}d_{i}$, we find:

(a1) Both (out)ingoing waves undergo a weak cascade, $\chi^{\pm}<1$, that
strengthens towards smaller scales. The trend is concurrent with tighter
scale-dependent dynamic alignment (SDDA) of fluctuations, a monotonic increase
in cross-helicity ($\sigma_{c}$), and a shift towards more negative residual
energy ($\sigma_{r}$) values

(a2) The ingoing waves transition to a strong cascade ($\chi^{-}\gtrsim 1$) at
$\lambda\approx 3\times 10^{4}d_{i}$; the associated spectral scalings deviate
from the expected weak-to-strong turbulence transitions. We explore the
possibility that “anomalous coherence’ effects may account for this
discrepancy in Section 6.

The domain canonically identified as the inertial range is comprised of two
distinct sub-inertial segments that exhibit distinct turbulence statistics.

For the subinertial range spanning $200-6000d_{i}$ and termed $R_{2}$ we find:

(b1) Spectral scaling indices for components parallel to the local mean field,
fluctuation (displacement), and perpendicular directions assume values of
$\alpha_{\ell_{||}}=-1.66\pm 0.05$, $\alpha_{\xi}=-1.56\pm 0.08$, and
$\alpha_{\lambda}=-1.49\pm 0.03$, respectively.

(b2) The “average” eddy assumes a field-aligned tube topology.

(b3) The alignment angle $\Theta^{ub}$ between velocity and magnetic-field
fluctuations monotonically increases towards smaller scales, while the
alignment $\Theta^{z}$ between the Elsässer fields remains roughly scale
independent ($\approx 35^{\circ}$). In both cases, an inverse relationship
between alignment angles and the intensity of field gradients is observed,
suggesting that the physical basis of alignment lies in the mutual shearing of
Elsässer fields during imbalanced collisions between counterpropagating wave
packets, as suggested in CSM15.

(b4) The cascade is strong for inwardly propagating waves ($\chi^{-}\gtrsim
1$) but weak for outwardly propagating ones, with $\chi^{+}$ increasing from
$0.1$ to $0.2$ as scales decrease from $\lambda\approx 10^{4}d_{i}$ to
$10^{2}d_{i}$.

(b5) The scaling exponents of the structure functions perpendicular to both
$\bm{B}_{\ell}$ and the fluctuation direction conform to the theoretical
models of CSM15 and MS17. However, the scaling profile in the parallel and
displacement components deviates from theoretical predictions, possibly due to
contamination from expansion effects.

For the subinertial range spanning $10-100d_{i}$ (termed $R_{1}$), we find:

(c1) The spectrum steepens, with spectral scaling indices for components
parallel to the local mean field, fluctuation (displacement), and
perpendicular directions assuming values of $\alpha_{\ell_{||}}=-1.97\pm
0.05$, $\alpha_{\xi}=-1.94\pm 0.06$, and $\alpha_{\lambda}=-1.64\pm 0.04$,
respectively.

(c2) A shift from isotropy in the plane perpendicular to $\bm{B}_{\ell}$
becomes evident, indicating a shift in eddy structures from tube-like to
ribbon-like, $\ell_{||}\gg\xi\gg\lambda$. While signatures of increasing SDDA
are observed, the result is potentially susceptible to errors in particle data
measurements.

(c3) The scaling exponents of the parallel and displacement components are a
linear function of order, while the perpendicular component exhibits a weakly
non-linear scaling profile. An overall transition towards “monofractal”
statistics and a weakening of intermittency, compared to $R_{2}$, are evident.

(d) At $\lambda\approx 8d_{i}$, the increase in aspect ratio ceases, and the
eddies transition to a quasi-isotropic state. This shift might be a signature
of the tearing instability, potentially leading to reconnection of the thin
current sheets, or it could result from turbulent energy being channeled into
an ion-cyclotron wave spectrum, consistent with the “helicity barrier” effect.

(e) The 2-point structure function method $SF^{n}_{2}$ is inadequate for
capturing the essential properties of the turbulent cascade, at smaller
scales. To accurately characterize steeper power laws at smaller spatial
scales, the use of a more sophisticated method such as the 5-point structure
function $SF^{n}_{5}$ is essential.

While our study doesn’t delve into the direct application of diagnostics for
expansion and imbalance effects, it’s interesting to note that preliminary
findings using data from the Wind mission show a notable correspondence with
the results from EBM simulations conducted by Shi et al. (2023). Specifically,
the scaling exponents in $R_{2}$ are consistent with CSM15 when
$\sigma_{c}\approx 1$, and become a linear function of order as the imbalance
decreases to $\sigma_{c}\approx 0$. Furthermore, extending the work of Verdini
et al. (2018) to higher order moments it is found that when intervals are
selected in such a way as to minimize the expansion effects the scaling
exponents in all three components are in striking agreement with those
predicted by the MS17 model (Sioulas et al. 2024, in progress).

In summary, our findings suggest that the models proposed by CSM15 and MS17,
which integrate SDDA as an intermittency effect and account for local
imbalance, possess the essential elements for a successful phenomenological
representation of imbalanced MHD turbulence. This assertion stems from several
key observations: firstly, the models provide scaling predictions for higher-
order moments in the perpendicular component of the magnetic field that align
well with our in-situ observations. Secondly, an inverse relationship between
alignment angles and the intensity of field gradients suggests that the
alignment mechanism originates from the mutual shearing of fields during
imbalanced collisions of wavepackets. However, it’s worth noting that certain
aspects of solar wind turbulence, such as the presence of two sub-inertial
ranges and anisotropic signatures, remain unaddressed by the models. This
suggests that incorporating additional effects, such as accounting for
inhomogeneity or the spherical polarization of fluctuations, could enhance the
models’ scaling predictions.

N.S. acknowledges useful conversations with Michael Stevens about the use of
SPC data. We acknowledge PSP/FIELDS team (PI: Stuart D. Bale, UC Berkeley) and
PSP/SWEAP team (PI: Justin Kasper, BWX Technologies) for the use of data. This
research was funded in part by the FIELDS experiment on the Parker Solar Probe
spacecraft, designed and developed under NASA contract NNN06AA01C; the NASA
Parker Solar Probe Observatory Scientist grant NNX15AF34G and the HERMES DRIVE
NASA Science Center grant No. 80NSSC20K0604. The instruments of PSP were
designed and developed under NASA contract NNN06AA01C. BC acknowledges the
support of NASA grant 80NSSC24K0171.

## References

  * Alexandrova et al. (2008) Alexandrova, O., Carbone, V., Veltri, P., & Sorriso-Valvo, L. 2008, ApJ, 674, 1153, doi: 10.1086/524056
  * Alfvén (1942) Alfvén, H. 1942, Nature, 150, 405, doi: 10.1038/150405d0
  * Andrés et al. (2019) Andrés, N., Sahraoui, F., Galtier, S., et al. 2019, Phys. Rev. Lett., 123, 245101, doi: 10.1103/PhysRevLett.123.245101
  * Angelopoulos et al. (2019) Angelopoulos, V., Cruce, P., Drozdov, A., et al. 2019, Space Science Reviews, 215, 9, doi: 10.1007/s11214-018-0576-4
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, ßr, 204, 49, doi: 10.1007/s11214-016-0244-5
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Nature, 1, doi: 10.1038/s41586-019-1818-7
  * Barnes & Hollweg (1974) Barnes, A., & Hollweg, J. V. 1974, Journal of Geophysical Research (1896-1977), 79, 2302, doi: https://doi.org/10.1029/JA079i016p02302
  * Batchelor et al. (1949) Batchelor, G. K., Townsend, A. A., & Jeffreys, H. 1949, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 199, 238, doi: 10.1098/rspa.1949.0136
  * Beresnyak (2011) Beresnyak, A. 2011, Phys. Rev. Lett., 106, 075001, doi: 10.1103/PhysRevLett.106.075001
  * Beresnyak (2012) Beresnyak, A. 2012, MNRAS, 422, 3495, doi: 10.1111/j.1365-2966.2012.20859.x
  * Beresnyak (2014) —. 2014, ApJ, 784, L20, doi: 10.1088/2041-8205/784/2/L20
  * Beresnyak (2015) —. 2015, ApJ, 801, L9, doi: 10.1088/2041-8205/801/1/L9
  * Beresnyak & Lazarian (2006) Beresnyak, A., & Lazarian, A. 2006, ApJ, 640, L175, doi: 10.1086/503708
  * Beresnyak & Lazarian (2008) —. 2008, ApJ, 682, 1070, doi: 10.1086/589428
  * Beresnyak & Lazarian (2009) —. 2009, ApJ, 702, 1190, doi: 10.1088/0004-637X/702/2/1190
  * Biskamp (2003) Biskamp, D. 2003, Magnetohydrodynamic Turbulence
  * Boldyrev (2002) Boldyrev, S. 2002, ApJ, 569, 841, doi: 10.1086/339403
  * Boldyrev (2005) Boldyrev, S. 2005, The Astrophysical Journal, 626, L37, doi: 10.1086/431649
  * Boldyrev (2006) —. 2006, Phys. Rev. Lett., 96, 115002, doi: 10.1103/PhysRevLett.96.115002
  * Boldyrev & Loureiro (2017) Boldyrev, S., & Loureiro, N. F. 2017, The Astrophysical Journal, 844, 125, doi: 10.3847/1538-4357/aa7d02
  * Boldyrev & Perez (2012) Boldyrev, S., & Perez, J. C. 2012, The Astrophysical Journal Letters, 758, L44, doi: 10.1088/2041-8205/758/2/L44
  * Bourouaine et al. (2020) Bourouaine, S., Perez, J. C., Klein, K. G., et al. 2020, ApJ, 904, L30, doi: 10.3847/2041-8213/abbd4a
  * Bowen et al. (2018) Bowen, T. A., Mallet, A., Bonnell, J. W., & Bale, S. D. 2018, ApJ, 865, 45, doi: 10.3847/1538-4357/aad95b
  * Bowen et al. (2020) Bowen, T. A., Bale, S. D., Bonnell, J. W., et al. 2020, Journal of Geophysical Research: Space Physics, 125, e2020JA027813, doi: https://doi.org/10.1029/2020JA027813
  * Bowen et al. (2020a) Bowen, T. A., Mallet, A., Bale, S. D., et al. 2020a, Phys. Rev. Lett., 125, 025102, doi: 10.1103/PhysRevLett.125.025102
  * Bowen et al. (2020b) Bowen, T. A., Mallet, A., Huang, J., et al. 2020b, ApJS, 246, 66, doi: 10.3847/1538-4365/ab6c65
  * Bowen et al. (2021) Bowen, T. A., Badman, S. T., Bale, S. D., et al. 2021, Nonlinear Interactions in Spherically Polarized Alfvénic Turbulence.  https://arxiv.org/abs/2110.11454
  * Bowen et al. (2023) Bowen, T. A., Bale, S. D., Chandran, B. D. G., et al. 2023, Mediation of Collisionless Turbulent Dissipation Through Cyclotron Resonance.  https://arxiv.org/abs/2306.04881
  * Bruno & Carbone (2013) Bruno, R., & Carbone, V. 2013, Living Reviews in Solar Physics, 10, 2, doi: 10.12942/lrsp-2013-2
  * Bruno et al. (2003) Bruno, R., Carbone, V., Sorriso-Valvo, L., & Bavassano, B. 2003, Journal of Geophysical Research (Space Physics), 108, 1130, doi: 10.1029/2002JA009615
  * Bruno et al. (2001) Bruno, R., Carbone, V., Veltri, P., Pietropaolo, E., & Bavassano, B. 2001, Planet. Space Sci., 49, 1201, doi: 10.1016/S0032-0633(01)00061-7
  * Burlaga (1991) Burlaga, L. F. 1991, J. Geophys. Res., 96, 5847, doi: 10.1029/91JA00087
  * Butterworth (1930) Butterworth, S. 1930, From 1923 to 1930, the journal was called Experimental Wireless and the Radio Engineer, 7, 536.  https://cir.nii.ac.jp/crid/1572261550684507392
  * Cerri & Califano (2017) Cerri, S. S., & Califano, F. 2017, New Journal of Physics, 19, 025007, doi: 10.1088/1367-2630/aa5c4a
  * Cerri et al. (2019) Cerri, S. S., Grošelj, D., & Franci, L. 2019, Frontiers in Astronomy and Space Sciences, 6, 64, doi: 10.3389/fspas.2019.00064
  * Cerri et al. (2022) Cerri, S. S., Passot, T., Laveder, D., Sulem, P. L., & Kunz, M. W. 2022, ApJ, 939, 36, doi: 10.3847/1538-4357/ac93fe
  * Chandran (2005) Chandran, B. D. G. 2005, Phys. Rev. Lett., 95, 265004, doi: 10.1103/PhysRevLett.95.265004
  * Chandran (2008) —. 2008, ApJ, 685, 646, doi: 10.1086/589432
  * Chandran (2018) Chandran, B. D. G. 2018, Journal of Plasma Physics, 84, 905840106, doi: 10.1017/S0022377818000016
  * Chandran & Perez (2019) Chandran, B. D. G., & Perez, J. C. 2019, Journal of Plasma Physics, 85, 905850409, doi: 10.1017/S0022377819000540
  * Chandran et al. (2015) Chandran, B. D. G., Schekochihin, A. A., & Mallet, A. 2015, ApJ, 807, 39, doi: 10.1088/0004-637X/807/1/39
  * Chandran et al. (2015) Chandran, B. D. G., Schekochihin, A. A., & Mallet, A. 2015, The Astrophysical Journal, 807, 39, doi: 10.1088/0004-637X/807/1/39
  * Chen (2016) Chen, C. H. K. 2016, Journal of Plasma Physics, 82, 535820602, doi: 10.1017/S0022377816001124
  * Chen et al. (2013) Chen, C. H. K., Bale, S. D., Salem, C. S., & Maruca, B. A. 2013, The Astrophysical Journal, 770, 125, doi: 10.1088/0004-637X/770/2/125
  * Chen et al. (2012) Chen, C. H. K., Mallet, A., Schekochihin, A. A., et al. 2012, ApJ, 758, 120, doi: 10.1088/0004-637X/758/2/120
  * Chen et al. (2011) Chen, C. H. K., Mallet, A., Yousef, T. A., Schekochihin, A. A., & Horbury, T. S. 2011, MNRAS, 415, 3219, doi: 10.1111/j.1365-2966.2011.18933.x
  * Chen et al. (2014) Chen, C. H. K., Sorriso-Valvo, L., Šafránková, J., & Němeček, Z. 2014, ApJ, 789, L8, doi: 10.1088/2041-8205/789/1/L8
  * Chhiber et al. (2021) Chhiber, R., Matthaeus, W. H., Bowen, T. A., & Bale, S. D. 2021, ApJ, 911, L7, doi: 10.3847/2041-8213/abf04e
  * Chhiber et al. (2020) Chhiber, R., Matthaeus, W. H., Oughton, S., & Parashar, T. N. 2020, Physics of Plasmas, 27, 062308, doi: 10.1063/5.0005109
  * Chhiber et al. (2020) Chhiber, R., Goldstein, M. L., Maruca, B. A., et al. 2020, The Astrophysical Journal Supplement Series, 246, 31, doi: 10.3847/1538-4365/ab53d2
  * Cho (2019) Cho, J. 2019, ApJ, 874, 75, doi: 10.3847/1538-4357/ab06f3
  * Cho & Lazarian (2003) Cho, J., & Lazarian, A. 2003, MNRAS, 345, 325, doi: 10.1046/j.1365-8711.2003.06941.x
  * Cho & Vishniac (2000) Cho, J., & Vishniac, E. T. 2000, ApJ, 539, 273, doi: 10.1086/309213
  * Comisso et al. (2018) Comisso, L., Huang, Y. M., Lingam, M., Hirvijoki, E., & Bhattacharjee, A. 2018, ApJ, 854, 103, doi: 10.3847/1538-4357/aaac83
  * Coppi et al. (1976) Coppi, B., Galvao, R., Pellat, R., Rosenbluth, M., & Rutherford, P. 1976, Soviet Journal of Plasma Physics, 2, 533
  * Cuesta et al. (2022) Cuesta, M. E., Chhiber, R., Roy, S., et al. 2022, ApJ, 932, L11, doi: 10.3847/2041-8213/ac73fd
  * D’Amicis et al. (2021) D’Amicis, R., Alielden, K., Perrone, D., et al. 2021, A&A, 654, A111, doi: 10.1051/0004-6361/202140600
  * D’Amicis et al. (2022) D’Amicis, R., Perrone, D., Velli, M., et al. 2022, Universe, 8, 352, doi: 10.3390/universe8070352
  * Davis et al. (2023) Davis, N., Chandran, B. D. G., Bowen, T. A., et al. 2023, ApJ, 950, 154, doi: 10.3847/1538-4357/acd177
  * Dmitruk et al. (2005) Dmitruk, P., Matthaeus, W. H., & Oughton, S. 2005, Physics of Plasmas, 12, 112304, doi: 10.1063/1.2128573
  * Dobrowolny et al. (1980) Dobrowolny, M., Mangeney, A., & Veltri, P. 1980, Phys. Rev. Lett., 45, 144, doi: 10.1103/PhysRevLett.45.144
  * Dong et al. (2022) Dong, C., Wang, L., Huang, Y.-M., et al. 2022, Science Advances, 8, eabn7627, doi: 10.1126/sciadv.abn7627
  * Dong et al. (2014) Dong, Y., Verdini, A., & Grappin, R. 2014, The Astrophysical Journal, 793, 118, doi: 10.1088/0004-637X/793/2/118
  * Duan et al. (2021) Duan, D., He, J., Bowen, T. A., et al. 2021, The Astrophysical Journal Letters, 915, L8, doi: 10.3847/2041-8213/ac07ac
  * Dudok de Wit et al. (2013) Dudok de Wit, T., Alexandrova, O., Furno, I., Sorriso-Valvo, L., & Zimbardo, G. 2013, Space Sci. Rev., 178, 665, doi: 10.1007/s11214-013-9974-9
  * Elsasser (1950) Elsasser, W. M. 1950, Phys. Rev., 79, 183, doi: 10.1103/PhysRev.79.183
  * Feller (1968) Feller, W. 1968, An Introduction to Probability Theory and Its Applications, Vol. 1 (Wiley).  http://www.amazon.ca/exec/obidos/redirect?tag=citeulike04-20{&}path=ASIN/0471257087
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, ßr, 204, 7, doi: 10.1007/s11214-015-0211-6
  * Franci et al. (2017) Franci, L., Cerri, S. S., Califano, F., et al. 2017, ApJ, 850, L16, doi: 10.3847/2041-8213/aa93fb
  * Franz & Harper (2002) Franz, M., & Harper, D. 2002, Planetary and Space Science, 50, 217, doi: https://doi.org/10.1016/S0032-0633(01)00119-2
  * Frisch (1995) Frisch, U. 1995, Turbulence: The Legacy of A.N. Kolmogorov, doi: 10.1017/CBO9781139170666
  * Frisch et al. (1978) Frisch, U., Sulem, P. L., & Nelkin, M. 1978, Journal of Fluid Mechanics, 87, 719, doi: 10.1017/S0022112078001846
  * Furth et al. (1963) Furth, H. P., Killeen, J., & Rosenbluth, M. N. 1963, The Physics of Fluids, 6, 459, doi: 10.1063/1.1706761
  * Galtier et al. (2000) Galtier, S., Nazarenko, S. V., Newell, A. C., & Pouquet, A. 2000, Journal of Plasma Physics, 63, 447, doi: 10.1017/S0022377899008284
  * Gerick et al. (2017) Gerick, F., Saur, J., & von Papen, M. 2017, The Astrophysical Journal, 843, 5, doi: 10.3847/1538-4357/aa767c
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, \apj, 438, 763, doi: 10.1086/175121
  * Goldstein et al. (1989) Goldstein, M. L., Roberts, D. A., & Matthaeus, W. H. 1989, Geophysical Monograph Series, 54, 113, doi: 10.1029/GM054p0113
  * Grappin (1986) Grappin, R. 1986, Physics of Fluids, 29, 2433, doi: 10.1063/1.865536
  * Grappin & Velli (1996) Grappin, R., & Velli, M. 1996, J. Geophys. Res., 101, 425, doi: 10.1029/95JA02147
  * Grappin et al. (1993) Grappin, R., Velli, M., & Mangeney, A. 1993, Phys. Rev. Lett., 70, 2190, doi: 10.1103/PhysRevLett.70.2190
  * Grappin et al. (2022) Grappin, R., Verdini, A., & Müller, W.-C. 2022, The Astrophysical Journal, 933, 246
  * Grauer et al. (1994) Grauer, R., Krug, J. H. A., & Marliani, C. 1994, Physics Letters A, 195, 335
  * Greco et al. (2018) Greco, A., Matthaeus, W. H., Perri, S., et al. 2018, ßr, 214, 1, doi: 10.1007/s11214-017-0435-8
  * Gurland & Tripathi (1971) Gurland, J., & Tripathi, R. C. 1971, The American Statistician, 25, 30.  http://www.jstor.org/stable/2682923
  * Higdon (1984) Higdon, J. C. 1984, ApJ, 285, 109, doi: 10.1086/162481
  * Hnat et al. (2011) Hnat, B., Chapman, S. C., Gogoberidze, G., & Wicks, R. T. 2011, Phys. Rev. E, 84, 065401, doi: 10.1103/PhysRevE.84.065401
  * Hollweg & Isenberg (2007) Hollweg, J. V., & Isenberg, P. A. 2007, Journal of Geophysical Research: Space Physics, 112, doi: https://doi.org/10.1029/2007JA012253
  * Horbury & Balogh (1997a) Horbury, T. S., & Balogh, A. 1997a, Nonlinear Processes in Geophysics, 4, 185, doi: 10.5194/npg-4-185-1997
  * Horbury & Balogh (1997b) —. 1997b, Nonlinear Processes in Geophysics, 4, 185, doi: 10.5194/npg-4-185-1997
  * Horbury & Balogh (2001) —. 2001, J. Geophys. Res., 106, 15929, doi: 10.1029/2000JA000108
  * Horbury et al. (2008) Horbury, T. S., Forman, M., & Oughton, S. 2008, \prl, 101, 175005, doi: 10.1103/PhysRevLett.101.175005
  * Horbury et al. (2018) Horbury, T. S., Matteini, L., & Stansby, D. 2018, MNRAS, 478, 1980, doi: 10.1093/mnras/sty953
  * Horbury et al. (2023) Horbury, T. S., Bale, S. D., McManus, M. D., et al. 2023, Physics of Plasmas, 30, 082905, doi: 10.1063/5.0123250
  * Howes (2015) Howes, G. G. 2015, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 373, 20140145, doi: 10.1098/rsta.2014.0145
  * Howes et al. (2012) Howes, G. G., Bale, S. D., Klein, K. G., et al. 2012, \apjl, 753, L19, doi: 10.1088/2041-8205/753/1/L19
  * Howes et al. (2008) Howes, G. G., Cowley, S. C., Dorland, W., et al. 2008, Journal of Geophysical Research (Space Physics), 113, A05103, doi: 10.1029/2007JA012665
  * Huang et al. (2023) Huang, Z., Sioulas, N., Shi, C., et al. 2023, ApJ, 950, L8, doi: 10.3847/2041-8213/acd7f2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iroshnikov (1963) Iroshnikov, P. S. 1963, Astronomicheskii Zhurnal, 40, 742.  https://ui.adsabs.harvard.edu/abs/1963AZh....40..742I
  * Kadomtsev & Pogutse (1967) Kadomtsev, B. B., & Pogutse, O. P. 1967, Soviet Journal of Experimental and Theoretical Physics, 24, 1172
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, ßr, 204, 131, doi: 10.1007/s11214-015-0206-3
  * Klein et al. (2012) Klein, K. G., Howes, G. G., TenBarge, J. M., et al. 2012, ApJ, 755, 159, doi: 10.1088/0004-637X/755/2/159
  * Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Kolmogorov (1962) Kolmogorov, A. N. 1962, Journal of Fluid Mechanics, 13, 82, doi: 10.1017/S0022112062000518
  * Kraichnan (1965) Kraichnan, R. H. 1965, The Physics of Fluids, 8, 1385
  * Landau & Lifshitz (1959) Landau, L. D., & Lifshitz, E. M. 1959, Fluid mechanics
  * Lithwick et al. (2007) Lithwick, Y., Goldreich, P., & Sridhar, S. 2007, ApJ, 655, 269, doi: 10.1086/509884
  * Loureiro & Boldyrev (2017) Loureiro, N. F., & Boldyrev, S. 2017, Phys. Rev. Lett., 118, 245101, doi: 10.1103/PhysRevLett.118.245101
  * Malagoli et al. (1996) Malagoli, A., Bodo, G., & Rosner, R. 1996, ApJ, 456, 708, doi: 10.1086/176691
  * Mallet & Schekochihin (2017) Mallet, A., & Schekochihin, A. A. 2017, MNRAS, 466, 3918, doi: 10.1093/mnras/stw3251
  * Mallet et al. (2015) Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2015, MNRAS, 449, L77, doi: 10.1093/mnrasl/slv021
  * Mallet et al. (2017) Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2017, Journal of Plasma Physics, 83, 905830609, doi: 10.1017/S0022377817000812
  * Mallet et al. (2016) Mallet, A., Schekochihin, A. A., Chandran, B. D. G., et al. 2016, MNRAS, 459, 2130, doi: 10.1093/mnras/stw802
  * Mallet et al. (2021) Mallet, A., Squire, J., Chandran, B. D. G., Bowen, T., & Bale, S. D. 2021, ApJ, 918, 62, doi: 10.3847/1538-4357/ac0c12
  * Mangeney (2001) Mangeney, A. 2001, in ESA Special Publication, Vol. 492, Sheffield Space Plasma Meeting: Multipoint Measurements versus Theory, ed. B. Warmbein, 53
  * Maron & Goldreich (2001) Maron, J., & Goldreich, P. 2001, ApJ, 554, 1175, doi: 10.1086/321413
  * Mason et al. (2006) Mason, J., Cattaneo, F., & Boldyrev, S. 2006, Phys. Rev. Lett., 97, 255002, doi: 10.1103/PhysRevLett.97.255002
  * Matteini et al. (2014) Matteini, L., Horbury, T. S., Neugebauer, M., & Goldstein, B. E. 2014, Geophysical Research Letters, 41, 259, doi: https://doi.org/10.1002/2013GL058482
  * Matteini et al. (2024) Matteini, L., Tenerani, A., Landi, S., et al. 2024, Physics of Plasmas, 31, 032901, doi: 10.1063/5.0177754
  * Matthaeus et al. (2008) Matthaeus, W. H., Pouquet, A., Mininni, P. D., Dmitruk, P., & Breech, B. 2008, Phys. Rev. Lett., 100, 085003, doi: 10.1103/PhysRevLett.100.085003
  * Matthaeus et al. (2012) Matthaeus, W. H., Servidio, S., Dmitruk, P., et al. 2012, ApJ, 750, 103, doi: 10.1088/0004-637X/750/2/103
  * Matthaeus & Velli (2011) Matthaeus, W. H., & Velli, M. 2011, Space Science Reviews, 160, 145, doi: 10.1007/s11214-011-9793-9