# Evolution of Magnetohydrodynamic Turbulence in the Expanding Solar Wind: Residual Energy and Intermittency

Chen Shi, Nikos Sioulas, Zesen Huang, Marco Velli, Anna Tenerani, Victor Réville

► To cite this version:

Chen Shi, Nikos Sioulas, Zesen Huang, Marco Velli, Anna Tenerani, et al. Evolution of Magnetohydrodynamic Turbulence in the Expanding Solar Wind: Residual Energy and Intermittency. The Astrophysical Journal, 2025, 979 (2), pp.152. 10.3847/1538-4357/ad9f38. hal-05398149

HAL Id: hal-05398149

https://hal.science/hal-05398149v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Evolution of Magnetohydrodynamic Turbulence in the Expanding Solar Wind: Residual Energy and Intermittency

Chen Shi (时辰)<fnref n="1" />, Nikos Sioulas<fnref n="1" />,<fnref n="2" />, Zesen Huang (黄泽森)<fnref n="1" />, Marco Velli<fnref n="1" />, Anna Tenerani<fnref n="3" />, and Victor Réville<fnref n="4" />

Received 2024 September 8; revised 2024 December 8; accepted 2024 December 12; published 2025 January 23

## Abstract

We conduct 3D magnetohydrodynamic simulations of decaying turbulence in the context of the solar wind. To account for the spherical expansion of the solar wind, we implement the expanding box model. The initial turbulence comprises uncorrelated counterpropagating Alfvén waves and exhibits an isotropic power spectrum. Our findings reveal the consistent generation of negative residual energy whenever nonlinear interactions are present, independent of the normalized cross helicity $\sigma_c$ and compressibility. The spherical expansion facilitates this process. The resulting residual energy is primarily distributed in the perpendicular direction, with $S_2(\mathbf{b}) - S_2(\mathbf{u}) \propto l_\perp$ or equivalently $-E_r \propto k_\perp^{-2}$. Here $S_2(\mathbf{b})$ and $S_2(\mathbf{u})$ are second-order structure functions of magnetic field and velocity respectively. In most runs, $S_2(\mathbf{b})$ develops a scaling relation $S_2(\mathbf{b}) \propto l_\perp^{1/2} (E_b \propto k_\perp^{-3/2})$. In contrast, $S_2(\mathbf{u})$ is consistently shallower than $S_2(\mathbf{b})$, which aligns with in situ observations of the solar wind. We observe that the higher-order statistics of the turbulence, which act as a proxy for intermittency, depend on the initial $\sigma_c$ and are strongly affected by the expansion effect. Generally, the intermittency is more pronounced when the expansion effect is present. Finally, we find that in our simulations, although the negative residual energy and intermittency grow simultaneously as the turbulence evolves, the causal relation between them seems to be weak, possibly because they are generated on different scales.

*Unified Astronomy Thesaurus concepts: Interplanetary turbulence (830); Magnetohydrodynamics (1964); Magnetohydrodynamical simulations (1966); Solar wind (1534)*

## 1. Introduction

It has long been observed that solar wind is a highly turbulent plasma system with fluctuations on a wide range of scales (see the review by R. Bruno & V. Carbone 2013, and references therein). Studying the solar wind turbulence is of great importance because turbulence is an important power source for the heating and acceleration of solar wind (S. R. Cranmer et al. 2007, 2015; A. Verdini et al. 2009; R. Lionello et al. 2014; A. Van Ballegooijen & M. Asgari-Targhi 2016; M. Shoda et al. 2019; V. Réville et al. 2020; N. Magyar & V. Nakariakov 2021; J. Halekas et al. 2023; Y. J. Rivera et al. 2024).

In the last decades, significant progress has been made on observations, numerical simulations, and theories of the solar wind turbulence. Satellite observations reveal that the turbulence in fast solar wind is usually highly Alfvénic, dominated by outward propagating Alfvén waves (J. Belcher & L. Davis 1971), while in slow solar wind Alfvénicity of the turbulence is typically lower than in the fast wind, but can be quite high in certain intervals (R. D'Amicis & R. Bruno 2015; R. D'Amicis et al. 2019), especially in the nascent solar wind as observed by Parker Solar Probe (O. Panasenco et al. 2020; T. Parashar et al. 2020).

Since compressible fluctuations are typically small in the solar wind, with $\delta n/n \lesssim 0.2$ (C. Shi et al. 2021) where n is the average plasma density and $\delta n$ is the fluctuation amplitude of the density, the solar wind turbulence is treated as an

incompressible magnetohydrodynamic (MHD) system in most theoretical and modeling works, where two Elsässer variables $z^\pm = u \mp b$, which are linear combinations of the velocity $u$ and magnetic field $b$ (in Alfvén speed) and represent the two counterpropagating Alfvén wave populations, are analyzed. A number of phenomenological models have been developed for the incompressible MHD turbulence. The weak ($|\delta b|/B \ll 1$), isotropic (in $k$-space), balanced ($z^+ \sim z^-$) turbulence model (P. Iroshnikov 1964; R. H. Kraichnan 1965) predicts a 1D power spectrum $E_{1D} \propto k^{-3/2}$. The weak, anisotropic, balanced model (P. Goldreich & S. Sridhar 1997) predicts a 1D power spectrum $E_{1D} \propto k_{\perp}^{-2}$. For strong, anisotropic turbulence, “critical balance” theory (P. Goldreich & S. Sridhar 1995), which balances the timescale for linear propagation and that for nonlinear eddy turnover, predicts $E_{1D} \propto k_{\perp}^{-5/3}$. Based on the critical balance theory, the scale-dependent dynamic alignment model (S. Boldyrev 2005; J. C. Perez & S. Boldyrev 2007) allows the spectral slope to be variable depending on how much the two Elsässer variables are aligned with each other. Imbalanced turbulence is more difficult to describe than balanced turbulence, and no simple phenomenological model has been established so far (M. Dobrowolny et al. 1980a, 1980b; R. Grappin et al. 1983; Y. Lithwick & P. Goldreich 2003; Y. Lithwick et al. 2007; A. Beresnyak & A. Lazarian 2010). Incompressible simulations conducted by J. C. Perez & S. Boldyrev (2009) show that in strong turbulence the two imbalanced Elsässer variables may have similar power spectra in spite of different amplitudes, while simulations conducted by A. Beresnyak & A. Lazarian (2009) show that the two Elsässer variables have very different structures.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, Los Angeles, CA 90095, USA; cshi1993@ucla.edu</footnote>

<footnote n="2">²Space Sciences Laboratory, University of California, Berkeley, Berkeley, CA 94720, USA</footnote>

<footnote n="3">³Department of Physics, The University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote n="4">⁴IRAP, Université Toulouse III–Paul Sabatier, CNRS, CNES, Toulouse, France</footnote>
Although the phenomenological models and previous numerical simulations have successfully explained some of the satellite observations, many mysteries still remain. One of the most outstanding problems is the prevailing negative residual energy, i.e., an excess of magnetic energy over the kinetic energy, in the solar wind turbulence (C. Chen et al. 2013, 2020; C. Shi et al. 2021; N. Sioulas et al. 2023). Many theoretical works have been carried out (W.-C. Müller & R. Grappin 2005; N. Yokoi & F. Hamba 2007; Y. Wang et al. 2011; S. Boldyrev et al. 2012a; G. Gogoberidze et al. 2012; G. G. Howes & K. D. Nielson 2013; S. Dorfman et al. 2024) to explain the generation of negative residual energy but they are not fully self-consistent and do not give consistent results, e.g., on the spectral slope of the residual energy. In addition, although most of the phenomenological models assume self-similarity, intermittency plays a nonnegligible role in MHD turbulence as it undermines the assumption of self-similarity (B. D. Chandran et al. 2015; A. Mallet & A. A. Schekochihin 2017; H. Wu et al. 2023). Besides, intermittency is an important means of energy dissipation and is observed to be directly correlated with plasma heating in the solar wind (N. Sioulas et al. 2022b, 2022c; C. Phillips et al. 2023).

In this study, we investigate, through 3D MHD simulations, the evolution of turbulence in the context of the solar wind with a focus on the residual energy and intermittency. The expanding box model (EBM; R. Grappin & M. Velli 1996; Y. Dong et al. 2014; A. Tenerani & M. Velli 2017; C. Shi et al. 2020, 2022; R. Grappin et al. 2022) was implemented because the spherical expansion of the solar wind may significantly change the evolution of turbulence as it leads to anisotropic decay of different components of the magnetic field and velocity and may result in mode conversion between different wave modes (Z. Huang et al. 2022). The paper is organized as follows. In Section 2, we describe the simulation setup. In Section 3, we present the simulation results. In Section 4 we discuss the relation between the residual energy and intermittency. In Section 5, we summarize this work.

## 2. Simulation Setup

We use the LAPS code,<fnref n="5" /> which is a 3D pseudospectral compressible MHD code with the EBM, to conduct the simulations. The algorithm of the code is described in detail by C. Shi et al. (2024a). In all the simulations, the domain is a rectangular box with initial size ($5R_s$)$^3$ ($R_s$ is the solar radius) and grid number 512$^3$. Besides dealiasing in k-space, explicit resistivity $\eta = 2 \times 10^{-5}$ and viscosity $\nu = 2 \times 10^{-5}$ are implemented to maintain numerical stability. We note that, because the code is based on the MHD equation in conserved form, the viscosity is implemented as $\partial_t(\rho u)_k \sim -k^2\nu(\rho u)_k$ where $(\rho u)_k$ is Fourier mode k of the conserved variable $\rho u$.

The initial fields consist of a uniform background and fluctuations. The background fields are $\rho_0 = B = 1$, $P_0 = 0.1006$, with normalization units $\bar{n} = 200 \text{ cm}^{-3}$, $\bar{B} = 250 \text{ nT}$, and subsequently $\bar{P} = \bar{B}^2/\mu_0 = 49.7 \text{ nPa}$ where $\mu_0$ is the permeability. The background magnetic field is within the equatorial plane (x-y plane) and has an angle of 8.1 with respect to the radial direction ($\hat{e}_x$), so that in simulations with expansion this angle increases to about 45° at 1 au. The adiabatic index is $\gamma = 1.5$ instead of 5/3 to prevent the plasma temperature from cooling too fast in the runs with expansion.

This choice only slightly modifies the thermodynamics and is not expected to impact our results significantly.

Fluctuations of velocity and magnetic field with 3D isotropic power spectra are added on top of the background fields. These initial fluctuations are added on the largest scales contained in the simulation domain such that $|k| \in [1/L, 32/L]$ where $L = 5R_s$ is the domain size. The reduced 1D spectra of the initial fluctuations roughly follow $|k|^{-1.3}$. Consequently, we are not able to observe the shallow "1/f range" that is usually observed in the solar wind (e.g., L. Matteini et al. 2018; Z. Huang et al. 2023) and may be generated due to the inverse cascade (B. D. Chandran 2018; R. Meyrand et al. 2023). We note that there are no forcing terms in the model equations, i.e., the turbulence is decaying. One should be aware that, with the expansion effect, in addition to dissipation, the turbulence also decays due to energy exchange with the background solar wind. The initial fluctuations are Alfvénic: for any wave mode $k$, there is $b_k \propto k \times b$, where $b_k$ is the magnetic field fluctuation of wavevector $k$. Usually, we use the normalized cross helicity $\sigma_c$ and normalized residual energy $\sigma_r$ to measure the Alfvénicity of the turbulence, and they are defined as

$$ \sigma_c = \frac{E_+ - E_-}{E_+ + E_-}, \quad \sigma_r = \frac{E_k - E_b}{E_k + E_b} \quad (1) $$

where $E_\pm$ represent the energy of the outward/inward Alfvén waves ($z^\pm$) and $E_{k,b}$ represent the kinetic and magnetic energies of the fluctuations. We note that $\sigma_c$ measures the correlation between the velocity and magnetic field, and $\sigma_r$ measures the correlation between the two Elsässer variables. At initialization, we control $\sigma_c$ by varying the correlation between the velocity fluctuation and magnetic field fluctuation, and we keep $\sigma_r$ exactly zero.

The rms of the magnetic field fluctuation is $b_{rms}/B \approx 0.14$ for all the runs. Thus, the nonlinear eddy turnover time is estimated to be $\tau_{nl} \sim L/2\pi b_{rms} \approx 5.7$ and the effective Reynolds number $Re$ (and Lundquist number $S$) is $Re = S \approx Lb_{rms}/\nu \approx 3.5 \times 10^4$. We note that because the background plasma has quite low $\beta$ ($\approx 0.2$) as we want the configuration to be close to the realistic solar wind in the inner heliosphere (A. Artemyev et al. 2022), the fluctuation level cannot be too strong, otherwise the simulation will be unstable due to formation of shocks. In contrast, Y. Dong et al. (2014) and R. Grappin et al. (2022) added strong turbulence (with $b_{rms}/B \approx 1$) to their simulations by adopting large $\beta$. The initial turbulence Mach number in our simulations is $M_s = u_{rms}/C_s \approx 0.36$, where $u_{rms} = b_{rms} = 0.14$ and $C_s = \sqrt{\gamma P_0/\rho_0} = 0.39$ is the sound speed.

In the runs with expansion, the initial radial location of the simulation domain is $R_0 = 30R_s$, and the radial speed of the box is $U_r = 1.167$ with normalization unit $\bar{U} = \bar{B}/\sqrt{\mu_0 m_p \bar{n}} = 385.6 \text{ km s}^{-1}$ where $m_p$ is the proton mass. We carry out six compressible-MHD runs, which are divided into three groups: Runs (10E, 10NE), Runs (05E, 05NE), and Runs (00E, 00NE). Here "E" and "NE" stand for "expansion" and "no expansion" respectively. Runs 10 have $\sigma_{c,0} = 1$, Runs 05 have $\sigma_{c,0} \approx 0.5$, and Runs 00 have $\sigma_{c,0} \approx 0$ where $\sigma_{c,0}$ is the initial normalized cross helicity. In addition, based on Runs 05NE and Run 00NE, we carry out two extra runs using the incompressible version of the LAPS code. These two runs are labeled Run 05IC and Run 00IC.

<footnote n="5">⁵ https://github.com/chenshihelio/LAPS</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Time evolution of $|\delta\mathbf{B}|/|\mathbf{B}|$ (blue), $\delta|\mathbf{B}|/|\mathbf{B}|$ (orange), and $\delta\rho/\rho$ (green) in different runs. Here $\delta\mathbf{b}$ is the rms of the magnetic field vector, $\delta|\mathbf{B}|$ is the rms of the magnetic field magnitude, and $\delta\rho$ is the rms of the density. $|\mathbf{B}|$ is the amplitude of the average magnetic field, and $\rho$ is the average density. Left, middle, and right panels show Runs 10, Runs 05, and Runs 00 respectively. The black dashed-dotted lines show the growth of $|\delta\mathbf{B}|/|\mathbf{B}|$ predicted by WKB theory, i.e., $\propto R^{1/2}$.</figcaption></figure>

# 3. Result

## 3.1. Compressibility

Figure 1 shows the time evolution of the normalized fluctuation levels of the magnetic field vector (blue), magnetic field magnitude (orange), and density (green) in different runs. Here the fluctuation level is defined as the rms of a specific quantity. We use the magnitude of the average magnetic field and the average density for normalization. In each panel, the solid curves correspond to the run with expansion, the dashed curves to the run without expansion, and the dotted curves to the incompressible run.

Comparing the three panels, we find that the initial $\sigma_c$ does not have a strong impact on the evolution of the three parameters. In runs with expansion, $|\delta\mathbf{B}|/|\mathbf{B}|$ increases at the beginning and then slowly decays after saturation. The increase is attributed to the slower decay of the Alfvén wave amplitude than the background magnetic field in the expanding solar wind (J. Belcher 1971; J. V. Hollweg 1974). However, we note that the increase is slower than the WKB prediction, i.e., $\propto R^{1/2}$, as shown by the black dashed-dotted curves in Figure 1, because of the nonlinear energy cascade. The later decrease is because of the dissipation of the turbulence energy, as can be seen in the runs without expansion. Because of the initial pressure imbalance, compressible fluctuations are generated soon after the simulations start. Interestingly, the density fluctuation is stronger in runs with lower $\sigma_c$, indicating that the nonlinear interaction plays an important role in the generation of compressible fluctuations. Similar to $|\delta\mathbf{B}|/|\mathbf{B}|$, the normalized density fluctuation increases with time in runs with expansion but decreases in runs without expansion. The magnetic compressibility $\delta|\mathbf{B}|/|\mathbf{B}|$ follows a similar trend with $\delta\rho/\rho$, i.e., it increases with time in runs with expansion and decreases in runs without expansion. We note that, in a recent study by L. Matteini et al. (2024), it was found that the magnetic field magnitude evolves toward uniform in 2D EBM hybrid simulations of balanced turbulence. This, however, is not observed in our MHD simulations, which show that the magnetic compressibility continues to increase with time. This discrepancy may imply that kinetic physics is necessary to produce the spherically polarized Alfvén waves that dominate the solar wind turbulence. In the EBM-MHD simulations conducted by J. Squire et al. (2020), the magnetic compressibility is observed to be small as the fluctuation level of the

magnetic field vector becomes similar to the background magnetic field, accompanied by the generation of rotational discontinuities (B. J. Vasquez & J. V. Hollweg 1998). However, our simulations are not directly comparable to those by J. Squire et al. (2020) because the amplitude of the fluctuations is relatively small. A detailed analysis of the nature of the compressible fluctuations in the EBM-MHD simulations is necessary but will be left for a future study.

## 3.2. Evolution of $\sigma_c$ and $\sigma_r$

In Figure 2, we show the time evolution of $\sigma_c$ (top) and $\sigma_r$ (bottom) in Runs 10 (a), Runs 05 (b), and Runs 00 (c). Blue lines with circles are runs with expansion, orange lines with triangles are runs without expansion, and green lines with squares are incompressible runs without expansion. With expansion, $|\sigma_c|$ gradually decreases in imbalanced turbulence (Runs 10E and 05E), because of the reflection of the outward propagating Alfvén waves due to the inhomogeneity of the background fields (M. Heinemann & S. Olbert 1980; M. Velli et al. 1991). Without expansion, Run 10NE does not evolve because nonlinear interaction is absent from the exactly Alfvénic status ($\sigma_{c,0} = 1$). In Run 05NE and Run 05IC, $\sigma_c$ increases with time, possibly because of the “dynamic alignment” (M. Dobrowolny et al. 1980a, 1980b), i.e., an initially imbalanced turbulence tends to evolve toward purely Alfvénic status because the energy decay rates of the two counterpropagating Alfvén wave populations are similar. For the (nearly) balanced turbulence (Runs 00), whether $\sigma_c$ evolves toward positive or negative is very sensitive to the initial condition. In the simulations conducted here, the initial $\sigma_c$ is slightly negative. Hence, in Run 00NE and Run 00IC, $\sigma_c$ decreases to more and more negative values. In Run 00E, $\sigma_c$ remains negative but stays at very small absolute values. This is because of the competition between the dynamic alignment, which tends to increase $|\sigma_c|$, and the expansion effect, which tends to decrease $|\sigma_c|$. By comparing the green and orange lines, we can see that the evolution of $\sigma_c$ in the incompressible runs (Runs 05IC and 00IC) does not show big differences from that in the compressible runs (Runs 05NE and Runs 00NE).

$\sigma_r$ evolves toward negative values in all the runs except for Run 10NE, indicating that nonlinear interaction naturally generates negative residual energy (R. Grappin et al. 1983;
<figure><img src="image_5.png" /><figcaption>Figure 2. Top panels: evolution of $\sigma_c$ in Runs 10 (a1), Runs 05 (b1), and Runs 00 (c1). Blue circles are runs with expansion, orange triangles are runs without expansion, and green squares are incompressible-MHD runs without expansion. The bottom panels have the same format as the top panels but show the evolution of $\sigma_r$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. Evolution of kinetic energy (blue), magnetic energy (orange), and total energy (black) with high-resolution output for Run 05E (a), Run 05NE (b), and Run 05IC (c). The bottom row shows blow-ups of the yellow shaded regions in the top row. For Run 05E, we have multiplied the energies by $(R/R_0)^3$ to compensate for the energy decay due to expansion based on WKB theory.</figcaption></figure>

N. Yokoi & F. Hamba 2007; Y. Wang et al. 2011; S. Boldyrev et al. 2012a; G. Gogoberidze et al. 2012; G. G. Howes & K. D. Nielson 2013; S. Dorfman et al. 2024), consistent with the prevailing negative residual energy observed in the solar wind (C. Chen et al. 2013, 2020; C. Shi et al. 2021). In addition, it is clear that $\sigma_r$ decays faster and to more negative

values in runs with expansion than in runs without expansion. This may be attributed to the fact that the expansion-induced decay of magnetic energy is slower than that of kinetic energy (Y. Dong et al. 2014; C. Shi et al. 2022), i.e.,

$$b_{\perp}^2 \propto R^{-2}, \quad \rho u_{\perp}^2 \propto R^{-4} \qquad (2)$$
in the absence of any coupling between **b** and **u**. This effect is significant for nonpropagating, perpendicular modes with $k \perp B_0$ (R. Meyrand et al. 2023). Similar to $\sigma_c$, $\sigma_r$ in the incompressible runs does not differ much from the compressible runs.

As an example, in Figure 3, we show the time evolution of kinetic (blue) and magnetic (orange) energies in Run 05E, Run 05NE, and Run 05IC with high time resolution. Total energy is also shown as the black curves. For Run 05E we have multiplied the energies by $(R/R_0)^3$ to compensate for the energy decay due to solar wind expansion (the WKB theory, J. Belcher 1971). We observe high-frequency oscillations of $E_k$ and $E_b$, which are anticorrelated so that the total energy does not oscillate. This oscillation is clearly a result of the effect of wave propagation (Y. Wang et al. 2011). From panel (b) or (c), we can estimate the oscillation period is roughly $T \approx 2.5 \approx L_x/2B_0$, and from panel (a) we see that the period increases gradually because the expansion increases the crossing time of Alfvén waves through the simulation domain. Figure 3 indicates that the magnetic energy excess is built up over multiple wave crossing times.

### 3.3. Second-order Structure Functions

We then investigate the second-order structure functions of different fields. The qth-order structure function of a field $b(x)$ is defined as

$$S_q(b, l) = \langle |b(x+l) - b(x)|^q \rangle_x \quad (3)$$

where $l$ is a given spatial increment, and $\langle \rangle_x$ means ensemble average or equivalently average over the whole simulation domain. The second-order structure function $S_2$ measures the mean square value of the fluctuation on scale $l$. One important relation is that, if $S_2$ scales exponentially with the spatial increment such that $S_2 \propto l^\alpha$, the power spectrum of the field obeys the scaling $E \propto k^{-(\alpha+1)}$ (E. W. Montroll & M. F. Shlesinger 1982).

Figure 4 shows $S_2$ of velocity (blue), magnetic field (orange dashed), and negative residual energy (magnetic energy minus kinetic energy, green dotted) at the end of the simulations ($t=200$), when turbulence has evolved for a sufficient time. The left, middle, and right columns show structure functions with $l$ along $x$, $y$, and $z$. The rows from top to bottom are Runs 10E, 05E, 05NE, 05IC, 00E, 00NE, and 00IC. Run 10NE is not shown because of the absence of nonlinear evolution. Anisotropy among the three directions is clearly observed in all the runs. Along $l_x$ (radial and initially quasi-parallel), no clear power-law relation is established. Along $l_y$, an extended power-law part with a slope slightly smaller than 0.5 forms in all the runs for both **u** and **b**. Along $l_z$, we get $S_2(b) \propto l_z^{0.5}$ and $S_2(u)$ shallower than $S_2(b)$. This is consistent with satellite observations (C. Chen et al. 2020; C. Shi et al. 2021), though the spectral slopes from the simulations are systematically shallower than satellite observations, which show that the magnetic field and velocity power spectra have slopes of $-5/3$ and $-3/2$ on average (C. Chen et al. 2013). By comparing different rows, one can see that $S_2(b) \propto l_z^{0.5}$ holds for all the runs, independent of the initial $\sigma_c$ and the expansion effect, while $S_2(u)$ is clearly shallower in the runs with expansion than in runs without expansion. Residual energy is generated mainly along the $z$-direction, i.e., the direction perpendicular to both the background magnetic field and the radial direction. The

negative residual energy has a power-law scaling $S_2 \propto l_z$, i.e., $-E_r \propto k_z^{-2}$ in Runs 10E, 05E, 05NE, 05IC, and 00E. This spectral slope is consistent with the WIND observation (C. Chen et al. 2013) as well as the prediction given by the eddy-damped quasi-normal Markovian (EDQNM) model of isotropic MHD turbulence (R. Grappin et al. 1983),<fnref n="6" /> which is verified by 3D MHD simulations with zero mean magnetic field (R. Grappin et al. 2016) and is consistent with the model by S. Boldyrev et al. (2012a) for anisotropic strong balanced turbulence. However, we note that there is so far no self-consistent theory for the spectral slope of residual energy, and different models can give different results. For example, Y. Wang et al. (2011) show, through analytic calculation, that in the scenario of weak turbulence negative residual energy is produced and follows $-E_r \propto k_\perp^{-1}$. The EDQNM model for strong, anisotropic MHD turbulence (G. Gogoberidze et al. 2012) also shows that negative residual energy arises due to nonlinear interaction but its spectrum follows $-E_r \propto k_\perp^{-5/3}$. Our results reveal that the generation of residual energy is strongly anisotropic with the presence of a finite background magnetic field. Even with a small mean field, e.g., along the $y$-direction in Run 05NE, residual energy is barely observed, possibly due to the “Alfvén effect” that dissipates the residual energy (R. H. Kraichnan 1965; R. Grappin et al. 2016).

In Figure 5, we show $S_2(l)$ of $z^+$ (blue) and $z^-$ (orange dashed) at the end of different runs. Similar to Figure 4, for $l = l_\hat{e}_z$ (quasi-parallel to $B_0$), the structure function does not evolve much. For $l = l_\hat{e}_y$, the structure function is determined by a mix of parallel and perpendicular effects. Therefore, we will focus on the right column, i.e., $l = l_\hat{e}_z$. $S_2(z^+)$ has very similar shapes, i.e., with a slope slightly shallower than 0.5, in all the runs. The slope of $S_2(z^-)$, however, behaves very differently from that of $S_2(z^+)$. For Runs 05NE, 05IC, 00NE, and 00IC, the slope is roughly 0.5, but for runs with expansion, it is strongly affected by $\sigma_c$. $S_2(z^-) \propto l_z$ for Run 10E, $S_2(z^-) \propto l_z^{2/3}$ for Run 05E, and $S_2(z^-) \propto l_z^{0.5}$ for Run 00E. That is to say, the perpendicular spectrum of $z^-$ is steeper as the turbulence gets more imbalanced. In previous numerical works without the expansion effect (J. C. Perez & S. Boldyrev 2009; J. C. Perez et al. 2012), the spectral slope of $z^-$ is $E \sim k_\perp^{-1.5}$ for both balanced and imbalanced turbulence, consistent with our results, but the spectrum of $z^+$ is steeper than that of $z^-$. This inconsistency may be a result of the difference in the simulation setup, as J. C. Perez & S. Boldyrev (2009) and J. C. Perez et al. (2012) implement driving forces for the turbulence while our simulations contain decaying turbulence. In a recent MHD simulation of decaying strong turbulence (L. Yang et al. 2023), the spectral slopes of $z^+$ and $z^-$ are roughly $-5/3$ and $-1$ respectively. In addition, observation by Parker Solar Probe shows that the $z^+$ spectrum is mostly steeper than the $z^-$ spectrum (C. Shi et al. 2021). The discrepancy between our simulation result and these previous studies is still unclear. Nonetheless, we note that R. Grappin et al. (2022), through a comprehensive parametric study with EBM simulations, find that the spectral slopes of $z^+$ can be affected by various factors, including the initial spectral slopes and the strength of turbulence.

### 3.4. Higher-order Statistics of the Magnetic Field

It is well known that intermittency develops in MHD turbulence, generating local structures in magnetic field and

<footnote n="6">⁶ In R. Grappin et al. (1983), the sign of the residual energy is not defined.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Second-order structure functions of velocity $S_2(u)$ (blue) and magnetic field $S_2(b)$ (orange dashed), and the difference between them $S_2(b) - S_2(u)$ (green dotted) at the end of simulations ($t = 200$). Columns from left to right are $l = l_{e_x}$, $l = l_{e_y}$, and $l = l_{e_z}$. Rows from top to bottom are Runs 10E, 05E, 05NE, 05IC, 00E, 00NE, and 00IC.</figcaption></figure>

velocity. A useful measure of the intermittency is the scaling exponents of structure functions. For homogeneous turbulence without intermittency, distribution of the fluctuations is typically assumed to be self-similar across different spatial

scales, in which case the slope ("scaling exponent") of the qth-order structure function is a linear function of $q$, i.e., "monofractal." With intermittency, the slope becomes a nonlinear function of $q$, i.e., "multifractal," due to the change
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Second-order structure functions of Elsässer variables $S_2(z^+)$ (blue) and $S_2(z^-)$ (orange dashed) at the end of simulations ($t = 200$). Columns from left to right are $l = l_x$, $l = l_y$, and $l = l_z$. Rows from top to bottom are Runs 10E, 05E, 05NE, 05IC, 00E, 00NE, and 00IC.</figcaption></figure>

in the distribution of the fluctuations as we move toward small scales. Observations have revealed that the magnetic field fluctuations in the solar wind are typically multifractal (e.g., L. Sorriso-Valvo et al. 1999; J. C. Palacios et al. 2022; N. Sioulas et al. 2022a).

In Figure 6, we show the magnetic field structure functions $S_q(b, l)$ at the end of Run 10E as an example. For each curve, we apply linear fitting to the range $L/16 \le l \le L/2$, which is marked by the yellow shades, and get the scaling exponents. The result of linear fitting is shown by the black dashed lines.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Structure functions of magnetic field $S_q(b, l)$ for (a) $l = l_x \hat{e}_x$, (b) $l = l_y \hat{e}_y$, and (c) $l = l_z \hat{e}_z$ at the end of Run 10E ($t = 200$). Curves for $q = 1–6$ are plotted. Yellow shades mark the scale ranges used for fitting the structure functions. The result of linear fitting is shown by the black dashed lines.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Scaling exponents as functions of $q$ along $x$ (blue squares), $y$ (orange crosses), and $z$ (green triangles) directions at the end of different runs. In each panel, the red dashed line shows the initial status ($t = 0$). The three black lines are $q/3$, $q/4$, and $1 - q^{0.69}$ for reference.</figcaption></figure>

In Figure 7, we show the fitted scaling exponents at the end of different runs. As Runs 05IC and 00IC show similar results to Runs 05NE and 00NE, they are not shown in this figure. Here the blue curves with squares correspond to $l_x$, the orange curves with crosses to $l_y$, and the green curves with triangles to $l_z$. For reference, the black dashed line shows $q/3$, which is the Kolmogorov turbulence model, the black dotted line shows $q/4$, which is the Iroshnikov–Kraichnan turbulence model, and the black dashed–dotted line is $1 - q^{0.69}$, which is a multifractal intermittency model based on the assumption of

strong, balanced turbulence (B. D. Chandran et al. 2015). We note that the result for Run 10NE is unreliable because the structure functions barely develop a power-law in this run. Anisotropy among $x$, $y$, and $z$-axes as well as multifractality are clearly observed in all the runs. Inspecting the results for $l_z$, we find that the scaling exponents in Run 00NE roughly follow the prediction by B. D. Chandran et al. (2015), as expected. Surprisingly, in Run 10E (strongly imbalanced turbulence with expansion), the scaling exponents also follow the prediction by B. D. Chandran et al. (2015), which, however, is based on the
<figure><img src="image_11.png" /><figcaption>Figure 8. Time evolution of kurtosis(J) (blue circles), kurtosis(ω) (orange triangles), and σr (black squares) in different runs.</figcaption></figure>

assumption of balanced turbulence. Why Run 10E shows stronger multifractality than Runs 05E and 00E is still unclear. It implies a complex interplay between the effect of nonlinear interaction and the effect of expansion on the evolution of intermittency.

Besides the scaling exponents, another useful quantification of intermittency is the kurtosis of current density $J = |∇ × b|$ and vorticity $ω = |∇ × u|$, which measures the strength of the intermittent current sheets and vortices. The kurtosis of a variable, which quantifies the deviation of its probability distribution function from the Gaussian distribution, is defined as

$$ \text{kurtosis}(f) = \frac{\langle f^4 \rangle_x}{\langle f^2 \rangle_x^2} $$

where again $\langle·\rangle_x$ stands for average over the simulation domain. We calculate the two quantities at the end of each run and write them in Figure 7. In Figure 8, we show the time evolution of kurtosis(J) (blue circles) and kurtosis(ω) (orange triangles) in different runs. Again, Run 10NE does not show observable evolution due to the lack of nonlinear interaction. Runs 05IC and 00IC are not shown because they show quite similar results to Runs 05NE and 00NE. Comparing kurtosis in different runs, we find that intermittency is obviously stronger in runs with expansion, possibly because of the selective decay of different components of the magnetic field and velocity, which gives rise to small-scale structures (Y. Dong et al. 2014). The kurtosis grows much faster in Run 10E than in Runs 05E and 00E, and the growth of the kurtosis is slightly faster in Run 05E than in Run 00E. This clearly shows that the evolution of intermittency is affected by $σ_c$. Although the discrepancy between the evolution of kurtosis in Runs 05NE and 00NE is not very pronounced, we expect to observe a slower evolution if the imbalance ($σ_c$) continues to increase, as implied by the stationary kurtosis in Run 10NE (Figure 8(a2)). Hence, the dependence of the growth rate of kurtosis on $σ_c$ is different in the

expansion runs and the non-expansion runs. This is reminiscent of the result shown in Figure 7, that is, the multifractality is stronger in Run 10E than in Runs 05E and 00E, while it is stronger in Run 00NE than in Run 05NE. This phenomenon is not fully understood yet and may imply a complex competition between the expansion effect and the nonlinear interaction in generating/dissipating the intermittency. A theory of intermittency for imbalanced turbulence in the expanding solar wind and an observational study to compare the strength of intermittency in solar wind streams with different $σ_c$ will be necessary.

## 4. Discussion: Is Residual Energy Related to Intermittency?

From Figure 8, one immediately notices that kurtosis(ω) is smaller than kurtosis(J) in all the runs, implying magnetic structures are stronger than velocity structures, consistent with in situ measurements by WIND (T. A. Bowen et al. 2018) and previous reduced-MHD simulations (V. Zhdankin et al. 2016). Consequently, one may conjecture that the negative residual energy is related to the intermittency. T. A. Bowen et al. (2018) find that there is a negative correlation between kurtosis(J) and $σ_r$, which implies that the negative residual energy is likely related to the intermittent magnetic structures. In Figure 8, we use black squares to show the time evolution of $σ_r$. One can see that as the simulation goes, $σ_r$ decreases while the kurtosis of both J and ω increases.

Figure 8 suggests that intermittent structures and negative residual energy are simultaneously generated as the turbulence evolves. However, whether the negative residual energy is produced by these intermittent structures is still unclear. To verify this point, we divide the simulation domain at the last frame ($t = 200$) of each run evenly into 16 × 16 × 16 cubes, i.e., each cube contains 32 × 32 × 32 grid points. For each cube, we calculate kurtosis(J), kurtosis(ω), and $σ_r$. Left and right columns of Figure 9 show the probability distribution of the data points for Run 05E and Run 05NE respectively. Black
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Probability distribution of kurtosis(J)–σ<sub>r</sub> (top), kurtosis(ω)–σ<sub>r</sub> (middle), and kurtosis(ω)–kurtosis(J) (bottom) based on the last snapshot (t = 200) of Run 05E (left) and Run 05NE (right). The simulation domain is divided into 16³ equal-size cubes and these quantities are calculated for each individual cube. Black lines are the median values of the y-axis values against the binned x-axis values. For Run 05E, the Pearson correlation coefficients are C.C.(kurt(J)–σ<sub>r</sub>) = -0.03, C.C.(kurt(ω)–σ<sub>r</sub>) = -0.02, and C.C.(kurt(ω)–kurt(J)) = 0.51. For Run 05NE, the Pearson correlation coefficients are C.C.(kurt(J)–σ<sub>r</sub>) = 0.02, C.C.(kurt(ω)–σ<sub>r</sub>) = -0.06, and C.C.(kurt(ω)–kurt(J)) = 0.35.</figcaption></figure>

curves are the median values of the y-axis values against the binned x-axis values. The other runs show similar results and hence are not shown here. We calculate the Pearson correlation coefficient between each pair of parameters and these correlation coefficients are written in the figure caption. The bottom row of Figure 9 shows that there is a positive correlation between kurtosis(J) and kurtosis(ω), though the correlation coefficients (0.51 for Run 05E and 0.35 for Run 05NE) are not high. This indicates that the intermittent structures in magnetic field and velocity are often colocated, but not always. As shown by the top and middle rows, kurtosis (J) and kurtosis(ω) do not show significant correlation with σ<sub>r</sub> in either run, with correlation coefficients of nearly zero. Although the data points spread widely in σ<sub>r</sub>, the median values of the kurtosis are quite constant. That is to say, at intermediate scales (L/16), regions with negative residual energy do not necessarily correspond to intermittent structures. We have carried out the same analysis with the simulation domain divided into 8 × 8 × 8 cubes, i.e., for a larger spatial scale (~L/8), and the result (not shown here) is almost identical.

Thus, we conclude that, in our simulations, negative residual energy and intermittency are simultaneously generated as a result of the evolution of turbulence, but the causal relation between them seems to be weak. The reason is that the residual energy is concentrated at large scales (Figure 4) while the intermittent structures are generated at smaller scales. Our result indicates that the negative residual energy is more likely produced by the wave–wave interaction (S. Boldyrev et al. 2012b; G. G. Howes & K. D. Nielson 2013) combined with the expansion effect. Intermittent structures may slightly contribute to the negative residual energy, considering the discrepancy between kurtosis(ω) and kurtosis(J) in all the runs, but they are likely insignificant. Here we note that the fact that the velocity intermittency is weaker than the magnetic intermittency is potentially due to the asymmetry between the momentum equation and the induction equation (V. Zhdankin et al. 2016). One can easily show that, in a 2D incompressible MHD system, the nonlinear term z̄<sup>±</sup> · ∇z<sup>±</sup> tends to contribute oppositely to the curl of the two Elsässer variables ω<sup>±</sup> = ∇ × z<sup>±</sup>. Since J = 1/2(ω<sup>+</sup> − ω<sup>−</sup>) and ω = 1/2(ω<sup>+</sup> + ω<sup>−</sup>), this asymmetry may lead to stronger current density than vorticity.
However, a rigorous analysis of this problem is still lacking and will be left for future work.

One possible explanation of the paradox between our simulations and satellite observations by T. A. Bowen et al. (2018) is that T. A. Bowen et al. (2018) adopted large time windows (one hour) to calculate these parameters and thus mixed large and small scales. Nonetheless, we note that due to artificial effects such as limited spatial resolution and lack of kinetic physics, the MHD simulations cannot capture all processes happening in the real solar wind. Moreover, besides the statistical analysis presented here, it would be beneficial to conduct a careful case study of the intermittent structures and their influences on the turbulence properties in MHD simulations in the future.

## 5. Summary

We conducted a set of 3D MHD simulations of solar wind turbulence with intermediate strength ($|\delta b|/B \sim 0.14$). The initialized fluctuations consist of counterpropagating Alfvén waves and have zero residual energy and varying normalized cross helicity. The key results are summarized below.

1. Negative residual energy is always produced when nonlinear interaction takes effect, regardless of normalized cross helicity. The effect of spherical expansion facilitates the generation of negative residual energy.

2. The magnetic field and velocity spectra are anisotropic and evolve differently. The magnetic field spectrum has a quite universal perpendicular slope of $-3/2$ while the velocity spectrum is shallower. The negative residual energy is observed primarily in the perpendicular direction and has a spectrum $-E_r \propto k_\perp^{-2}$ in most runs.

3. The spectral slope (along the perpendicular direction) of $z^+$ (outward) is quite universal and slightly shallower than $-3/2$, while the spectral slope of $z^-$ (inward) depends strongly on $\sigma_c$ when the expansion effect is turned on such that the imbalanced turbulence has a steeper $z^-$ spectrum. Without expansion, the $z^-$ spectrum has a slope of $-3/2$ for both the balanced runs (Runs 00NE and 00IC) and imbalanced runs (Runs 05NE and 05IC).

4. Runs with the expansion effect generate stronger intermittent structures in both magnetic field and velocity than the runs without expansion. The evolution of intermittency depends on $\sigma_c$ but the correlation between the intermittency and $\sigma_c$ is different in runs with expansion and in runs without.

5. Growth of negative residual energy is accompanied by the generation of intermittent structures. However, the causal relation between the negative residual energy and intermittency seems to be weak.

We emphasize that the strength of turbulence in our simulations is smaller than what is observed in the young solar wind, where $|\delta b|/B$ often reaches unity and thus magnetic “switchbacks” may form (S. Bale et al. 2019; J. C. Kasper et al. 2019; A. Tenerani et al. 2020, 2021). In addition, in the solar wind, the fluctuations are typically spherically polarized with $|B| = \text{const}$ (L. Matteini et al. 2018, 2024). In numerical simulations, although there have been efforts to construct spherically polarized magnetic field in 3D (F. Valentini et al. 2019; Z. Johnston et al. 2022; J. Squire & A. Mallet 2022; C. Shi et al. 2024b), it is nontrivial to impose a constant-$|B|$ magnetic field with a specified spectral slope (D. A. Roberts 2012).

## Acknowledgments

This study is supported by NSF SHINE #2229566, NASA HTMS #80NSSC20K1275, and NASA ECIP #80 NSSC23K1064. The numerical simulations are conducted on Extreme Science and Engineering Discovery Environment (XSEDE) EXPANSE at San Diego Supercomputer Center (SDSC) through allocation No. TG-AST200031, which is supported by National Science Foundation grant number ACI-1548562 (J. Towns et al. 2014).

Software: Matplotlib (J. D. Hunter 2007), NumPy (C. R. Harris et al. 2020), LAPS (C. Shi et al. 2024a).

## ORCID iDs

Chen Shi (时辰) https://orcid.org/0000-0002-2582-7085

Nikos Sioulas https://orcid.org/0000-0002-1128-9685

Zesen Huang (黄泽森) https://orcid.org/0000-0001-9570-5975

Marco Velli https://orcid.org/0000-0002-2381-3106

Anna Tenerani https://orcid.org/0000-0003-2880-6084

Victor Réville https://orcid.org/0000-0002-2916-3837

## References

Artemyev, A., Shi, C., Lin, Y., et al. 2022, ApJ, 939, 85

Bale, S., Badman, S., Bonnell, J., et al. 2019, Natur, 576, 237

Belcher, J. 1971, ApJ, 168, 509

Belcher, J., & Davis, L. 1971, JGR, 76, 3534

Beresnyak, A., & Lazarian, A. 2009, ApJ, 702, 460

Beresnyak, A., & Lazarian, A. 2010, ApJL, 722, 110

Boldyrev, S. 2005, ApJL, 626, L37

Boldyrev, S., Perez, J. C., & Wang, Y. 2012a, in ASP Conf. Ser. 459, 6th Int. Conf., Numerical Modeling of Space Plasma Flows, ed. N. V. Pogorelov, J. A. Font, & G. P. Zank (San Francisco, CA: ASP), 3

Boldyrev, S., Perez, J. C., & Zhdankin, V. 2012b, in AIP Conf. Proc. 1436 (Melville, NY: AIP), 18

Bowen, T. A., Mallet, A., Bonnell, J. W., & Bale, S. D. 2018, ApJ, 865, 45

Bruno, R., & Carbone, V. 2013, LRSP, 10, 1

Chandran, B. D. 2018, JPIPh, 84, 905840106

Chandran, B. D., Schekochihin, A. A., & Mallet, A. 2015, ApJ, 807, 39

Chen, C., Bale, S., Bonnell, J., et al. 2020, ApJS, 246, 53

Chen, C., Bale, S., Salem, C., & Maruca, B. 2013, ApJ, 770, 125

Cranmer, S. R., Asgari-Targhi, M., Miralles, M. P., et al. 2015, RSPTA, 373, 20140148

Cranmer, S. R., Van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS, 171, 520

D'Amicis, R., & Bruno, R. 2015, ApJ, 805, 84

D'Amicis, R., Matteini, L., & Bruno, R. 2019, MNRAS, 483, 4665

Dobrowolny, M., Mangeney, A., & Veltri, P. 1980a, in Proc. of IAU Symp. 91, Solar and Interplanetary Dynamics, ed. M. Dryer & E. Tandberg-Hanssen (Dordrecht: D. Reidel Publishing), 143

Dobrowolny, M., Mangeney, A., & Veltri, P. 1980b, PhRvL, 45, 144

Dong, Y., Verdini, A., & Grappin, R. 2014, ApJ, 793, 118

Dorfman, S., Abler, M., Boldyrev, S., Chen, C., & Greess, S. 2024, arXiv:2409.20442

Gogoberidze, G., Chapman, S. C., & Hnat, B. 2012, PhPl, 19, 102310

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

Goldreich, P., & Sridhar, S. 1997, ApJ, 485, 680

Grappin, R., Leorat, J., & Pouquet, A. 1983, A&A, 126, 51

Grappin, R., Müller, W.-C., & Verdini, A. 2016, A&A, 589, A131

Grappin, R., & Velli, M. 1996, JGRA, 101, 425

Grappin, R., Verdini, A., & Müller, W.-C. 2022, ApJ, 933, 246

Halekas, J., Bale, S., Berthomier, M., et al. 2023, ApJ, 952, 26

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Heinemann, M., & Olbert, S. 1980, JGRA, 85, 1311

Hollweg, J. V. 1974, JGR, 79, 1539

Howes, G. G., & Nielson, K. D. 2013, PhPl, 20, 072302

Huang, Z., Shi, C., Sioulas, N., & Velli, M. 2022, ApJ, 935, 60

Huang, Z., Sioulas, N., Shi, C., et al. 2023, ApJL, 950, L8

Hunter, J. D. 2007, CSE, 9, 90

Iroshnikov, P. 1964, SvA, 7, 566

Johnston, Z., Squire, J., Mallet, A., & Meyrand, R. 2022, PhPl, 29, 072902

Kasper, J. C., Bale, S. D., Belcher, J. W., et al. 2019, Natur, 576, 228

Kraichnan, R. H. 1965, PhFl, 8, 1385
Lionello, R., Velli, M., Downs, C., et al. 2014, ApJ, 784, 120

Lithwick, Y., & Goldreich, P. 2003, ApJ, 582, 1220

Lithwick, Y., Goldreich, P., & Sridhar, S. 2007, ApJ, 655, 269

Magyar, N., & Nakariakov, V. 2021, ApJ, 907, 55

Mallet, A., & Schekochihin, A. A. 2017, MNRAS, 466, 3918

Matteini, L., Stansby, D., Horbury, T., & Chen, C. H. 2018, ApJL, 869, L32

Matteini, L., Tenerani, A., Landi, S., et al. 2024, PhPl, 31, 032901

Meyrand, R., Squire, J., Mallet, A., & Chandran, B. D. 2023, arXiv:2308.10389

Montroll, E. W., & Shlesinger, M. F. 1982, PNAS, 79, 3380

Müller, W.-C., & Grappin, R. 2005, PhRvL, 95, 114502

Palacios, J. C., Bourouaine, S., & Perez, J. C. 2022, ApJL, 940, L20

Panasenco, O., Velli, M., Damicis, R., et al. 2020, ApJS, 246, 54

Parashar, T., Goldstein, M., Maruca, B., et al. 2020, ApJS, 246, 58

Perez, J. C., & Boldyrev, S. 2007, ApJL, 672, L61

Perez, J. C., & Boldyrev, S. 2009, PhRvL, 102, 025003

Perez, J. C., Mason, J., Boldyrev, S., & Cattaneo, F. 2012, PhRvX, 2, 041005

Phillips, C., Bandyopadhyay, R., McComas, D. J., & Bale, S. D. 2023, MNRAS: Letters, 519, L1

Réville, V., Velli, M., Panasenco, O., et al. 2020, ApJS, 246, 24

Rivera, Y. J., Badman, S. T., Stevens, M. L., et al. 2024, Sci, 385, 962

Roberts, D. A. 2012, PhRvL, 109, 231102

Shi, C., Tenerani, A., Rappazzo, A. F., & Velli, M. 2024a, FrASS, 11, 1412905

Shi, C., Velli, M., Panasenco, O., et al. 2021, A&A, 650, A21

Shi, C., Velli, M., Tenerani, A., Rappazzo, F., & Réville, V. 2020, ApJ, 888, 68

Shi, C., Velli, M., Tenerani, A., Réville, V., & Rappazzo, F. 2022, ApJ, 928, 93

Shi, C., Velli, M., Toth, G., et al. 2024b, ApJL, 964, L28

Shoda, M., Suzuki, T. K., Asgari-Targhi, M., & Yokoyama, T. 2019, ApJL, 880, L2

Sioulas, N., Huang, Z., Shi, C., et al. 2023, ApJL, 943, L8

Sioulas, N., Huang, Z., Velli, M., et al. 2022a, ApJ, 934, 143

Sioulas, N., Shi, C., Huang, Z., & Velli, M. 2022b, ApJL, 935, L29

Sioulas, N., Velli, M., Chhiber, R., et al. 2022c, ApJ, 927, 140

Sorriso-Valvo, L., Carbone, V., Veltri, P., Consolini, G., & Bruno, R. 1999, GeoRL, 26, 1801

Squire, J., Chandran, B. D., & Meyrand, R. 2020, ApJL, 891, L2

Squire, J., & Mallet, A. 2022, JPIPh, 88, 175880503

Tenerani, A., Sioulas, N., Matteini, L., et al. 2021, ApJL, 919, L31

Tenerani, A., & Velli, M. 2017, ApJ, 843, 26

Tenerani, A., Velli, M., Matteini, L., et al. 2020, ApJS, 246, 32

Towns, J., Cockerill, T., Dahan, M., et al. 2014, CSE, 16, 62

Valentini, F., Malara, F., Sorriso-Valvo, L., Bruno, R., & Primavera, L. 2019, ApJL, 881, L5

Van Ballegooijen, A., & Asgari-Targhi, M. 2016, ApJ, 821, 106

Vasquez, B. J., & Hollweg, J. V. 1998, JGRA, 103, 335

Velli, M., Grappin, R., & Mangeney, A. 1991, GApFD, 62, 101

Verdini, A., Velli, M., Matthaeus, W. H., Oughton, S., & Dmitruk, P. 2009, ApJL, 708, 116

Wang, Y., Boldyrev, S., & Perez, J. C. 2011, ApJL, 740, L36

Wu, H., Huang, S., Wang, X., et al. 2023, ApJL, 947, L22

Yang, L., He, J., Verscharen, D., et al. 2023, NatCo, 14, 7955

Yokoi, N., & Hamba, F. 2007, PhPl, 14, 112904

Zhdankin, V., Boldyrev, S., & Uzdensky, D. A. 2016, PhPl, 23, 055705