# Long-lived Equilibria in Kinetic Astrophysical Plasma Turbulence

Mario Imbrogno, Claudio Meringolo, Sergio Servidio, Alejandro Cruz-Osorio,
Benoît Cerutti, Francesco Pegoraro

► To cite this version:

Mario Imbrogno, Claudio Meringolo, Sergio Servidio, Alejandro Cruz-Osorio, Benoît Cerutti, et al. "Long-lived Equilibria in Kinetic Astrophysical Plasma Turbulence. The Astrophysical Journal Letters, 2024, 972, 10.3847/2041-8213/ad6b9d. insu-04836809

HAL Id: insu-04836809

https://insu.hal.science/insu-04836809v1

Submitted on 15 Dec 2024

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>

# Long-lived Equilibria in Kinetic Astrophysical Plasma Turbulence

Mario Imbrogno¹, Claudio Meringolo², Sergio Servidio³, Alejandro Cruz-Osorio⁴, Benoît Cerutti⁵, and Francesco Pegoraro⁶,⁷

¹Dipartimento di Fisica, Università della Calabria, Arcavacata (Cosenza), 87036, Italy

²Institut für Theoretische Physik, Goethe Universität, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany

³Dipartimento di Fisica, Università della Calabria, Arcavacata di Rende, 87036, Italy

⁴Instituto de Astronomía, Universidad Nacional Autónoma de México, AP 70-264, 04510 Ciudad de México, Mexico

⁵Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

⁶Dipartimento di Fisica “Enrico Fermi”, Università di Pisa, Pisa, 56122, Italy

⁷INAF-IAPS, via Fosso del Cavaliere 100, Roma, 00133, Italy

Received 2024 May 24; revised 2024 August 2; accepted 2024 August 6; published 2024 August 22

## Abstract

Turbulence in classical fluids is characterized by persistent structures that emerge from the chaotic landscape. We investigate the analogous process in fully kinetic plasma turbulence by using high-resolution, direct numerical simulations in two spatial dimensions. We observe the formation of long-lived vortices with a profile typical of macroscopic, magnetically dominated force-free states. Inspired by the Harris pinch model for inhomogeneous equilibria, we describe these metastable solutions with a self-consistent kinetic model in a cylindrical coordinate system centered on a representative vortex, starting from an explicit form of the particle velocity distribution function. Such new equilibria can be simplified to a Gold-Hoyle solution of the modified force-free state. Turbulence is mediated by the long-lived structures, accompanied by transients in which such vortices merge and form self-similarly new metastable equilibria. This process can be relevant to the comprehension of various astrophysical phenomena, going from the formation of plasmons in the vicinity of massive compact objects to the emergence of coherent structures in the heliosphere.

*Unified Astronomy Thesaurus concepts: Space plasmas (1544); High energy astrophysics (739); Plasma astrophysics (1261)*

## 1. Introduction

Astrophysical turbulence remains among the most fascinating phenomena, characterizing diverse systems ranging from the heliosphere to interstellar medium and compact object environments (Goldstein et al. 1995; Baiotti & Rezzolla 2017; Ripperda et al. 2022). This process covers a wide range of length- and timescales, from large-scale eddies to subelectron scales (Sahraoui et al. 2009). In this scenario, the cascade process is envisioned as a flux of energy from large-scale shears and boundary layers to the scales typical of particle interactions, where energy conversion is taking place (Matthaeus et al. 2015).

Turbulence is generally synonymous with randomness and unpredictability, although this is not quite correct. Persistent, long-lived structures indeed emerge from such a chaotic state, as observed in several systems (Chavanis & Sommeria 1998; Alexandrova 2008; Karimabadi et al. 2013). This “zoo” of coherent patterns can be qualitatively cataloged between vortices, waves, and discontinuous layers (Matthaeus et al. 2015). In this context, it is also important to mention large-scale structures such as Alfvén vortices (Pokhotelov & Petviashvili 1992; Alexandrova 2008).

Although considerable effort has been devoted to the characterization of persistent structures in classical (viscous) fluids (Montgomery et al. 1992), much less is known about the collisionless, magnetized counterpart. In plasma turbulence, stable structures may permeate the system and travel

undisturbed over long timescales, as one would expect in stellar winds (Borovsky 2008; Pecora et al. 2019) and accretion flows (Ripperda et al. 2020; Nathanail et al. 2022). These patterns, known as “plasmons” or “magnetic vortices” for their geometrical resemblance to hydrodynamical vortices, might be a crucial element of particle energization and dissipation (Drake et al. 2010; Petropoulou et al. 2016; Khabarova et al. 2021; El Mellah et al. 2022). Despite intense investigations, very little is known about their internal structure (Allanson et al. 2016; Lukin et al. 2018), mostly because of the coupling between large scales and characteristic plasma length scales.

In this Letter, we describe the process of coherent structure formation in fully kinetic plasma turbulence by using numerical simulations in 2.5D (2D in space, with 3D field components). We observe the formation of long-lived coherent structures typical of macroscopic, magnetically dominated force-free states. These metastable solutions can be described with a self-consistent kinetic model in a cylindrical coordinate frame, starting from an explicit form of the particle velocity distribution function. Such new equilibria can be simplified to a modified force-free state whose description can have several applications in all those studies concerned with the formation of coherent structures in astrophysical plasmas, such as plasmons in accretion flows and persistent flux ropes in the solar wind.

## 2. Methods

Our simulations are based on a full kinetic model of relativistic plasma by using the well-tested particle-in-cell (PIC) code Zeltron (Cerutti et al. 2013), which solves the equations of motion for a distribution of charged particles (i.e.,

<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.
characteristic curves for the Vlasov equation) when coupled to Maxwell's equations expressed in terms of a total magnetic field **b**, an electric field **ε**, a current density **j**:=∑<sub>α</sub>q<sub>α</sub>n<sub>α</sub>**u**<sub>α</sub>, and a charge density ρ<sub>c</sub>:=∑<sub>α</sub>q<sub>α</sub>n<sub>α</sub>, with n<sub>α</sub> the number density of each species, **u**<sub>α</sub> the bulk velocity, q<sub>α</sub> the charge, and α the species index representing either protons (p) or electrons (e).

The simulation setup closely follows that presented in Meringolo et al. (2023), with $N_x = N_y = 16,384$ mesh points in a square of side $L_0 \approx 5461d_e$, $d_e$ being the electron skin depth. We employ a realistic mass ratio $m_p/m_e = 1836$ and a total number of $\approx 2.7 \times 10^9$ macroparticles. Specific details on the simulation and the parameters are given in Appendix A. We impose large-scale random initial conditions for the magnetic field with a superposition of low-wavenumber Fourier modes so as to achieve a strong turbulence state (Meringolo et al. 2023). In particular, we set $\delta b/b_{0z} \sim 1$, where $\delta b$ is the rms of the magnetic field fluctuations and $b_{0z}$ is the out-of-plane (along z) mean magnetic field strength. Given the size of the magnetization $\sigma \approx 1$, which expresses the ratio between the magnetic pressure and the enthalpy density (the latter is the sum of the mass-energy density $\sum_{\alpha} n_{\alpha} m_{\alpha} c^2$ and the internal energy density $\sum_{\alpha} \rho_{\alpha} (\Gamma_{\alpha} - 1)^{-1}$ weighted by the adiabatic index $\Gamma_{\alpha}$, where $P_{\alpha}$ is the partial pressure related to the $\alpha$th species), and the plasma beta $\beta_p = \beta_e = 3 \times 10^{-3}$, which quantifies the ratio between kinetic and magnetic pressure, the observed dynamics is that of a weakly relativistic plasma.

### 3. Turbulence and Coherent Structures

The initial conditions rapidly produce a turbulent cascade, similar to fluid plasma models (Servidio et al. 2009). We concentrate our analysis when the peak of the averaged current density, which also corresponds to the time of the most intense nonlinear activity (Servidio et al. 2012), has long been reached (see the end of Appendix A for further information). Furthermore, under these conditions, a balance between the large-scale energy flux and the collisionless energy-conversion mechanisms is established, thus yielding a quasi-steady state. The power spectrum at that time is fully developed, as reported in Figure 1(a), and is consistent with typical observations of astrophysical turbulence (Bale et al. 2005; Alexandrova et al. 2009; Sahraoui et al. 2009).

As the broadband turbulence develops, we observe the appearance of coherent structures that move through the turbulent background, as can be inferred from the out-of-plane component of the magnetic potential $a_z$ in Figure 1(b). Similarly to MHD (Matthaeus et al. 2008; Servidio et al. 2008), kinetic turbulence tends to form local correlations, where the current **j** manifests a substantial tendency to align with the total magnetic field **b**. We estimate its strength by computing the $\cos(j, b) = j \cdot b / |j| |b|$, and as displayed in Figure 1(b), dominant structures are visible, whose morphology resembles hydrodynamic swirls and cyclones. They have a clear “eye,” where the alignment is net, and an outside region with advecting arms, where the cosine regularly changes sign while maintaining an overall circular symmetry. The characteristic size of such structures is $\approx 140d_e(3d_p)$ for the eye and $\approx 515d_e(12d_p)$ for the spiral arms, typical of inertial range turbulence, as can be seen from Figure 1(a). Note that the largest vortices (eye and arms) are on the order of a few $d_p$, corresponding to the correlation length, $\lambda_c$, set by our initial

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

<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) Magnetic power spectrum at a time beyond the peak of nonlinear activity vs. $kd_e$. The vertical (orange) line represents the persistent structure wavenumber $k_0 d_e$, which identifies the typical scale of the vortices. (b) 2D contour of the magnetic vector potential component $a_z$ (right side) and the cosine angle between the current density and the total magnetic field (left side). The y-axis is normalized to both the electron (left side) and the proton skin depth $d_p$ (right side).</figcaption></figure>

conditions of homogeneous turbulence, which, for the analyzed configuration, is $\lambda_c \approx 440d_e \approx 10d_p$.

To gain insight into the properties of these long-lived features, we focus our attention on regions where the magnitude of the magnetic potential $a_z$ at the center of the structure is much larger than its own rms (Servidio et al. 2009). After the time of maximum turbulence, we select four main vortices, which are marked with (white) circumferences in Figure 1(b). For each of these, we consider a local cylindrical coordinate system centered at the O-point of the vector potential and use it to produce locally azimuthally averaged quantities. More precisely, for any generic field $h$, we first transform it to the cylindrical coordinate frame, i.e., $h(x, y) \rightarrow h(r, \phi)$, and then compute the corresponding average as $H(r):=(2\pi)^{-1} \int_{0}^{2\pi} h(r, \phi') d\phi'$ (hereafter, we use capital letters to indicate azimuthally averaged quantities).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Radial behavior of the angle between the current density and the total magnetic field (azimuthal averages) for each structure depicted in Figure 1(b). (b) Force-free parameter $\lambda(r)$. We define $r_0$ as the vortex eye radius, namely, the position of the first spiraling arm.</figcaption></figure>

In Figure 2(a), we show the angle between the azimuthally averaged current density **J** and the magnetic field **B**, namely, $\Theta = \cos^{-1}(\|\mathbf{J} \cdot \mathbf{B}/|\mathbf{J}|\|\mathbf{B}|)$, evaluated inside the structures. The alignment here is extremely high, and all the vortices, which vary in size and intensity, exhibit the same behavior: they form highly aligned cores (eyes) and wall boundaries followed by spiraling arms (McWilliams 1984; Carnevale et al. 1991; Powell & Houston 1998; Haller 2005; Servidio et al. 2010a). Such spirals manifest alternate **J**–**B** alignment, implying a characteristic radial mode. We denote $r_0$ as the radius of the vortex eye, defined as the distance from the center to the first spiraling arm where the scalar product **J**·**B** changes sign (see Figure 2). The alignment is progressively weakened, and the plasmoid then blends with the background on scales $O(4r_0)$, which qualitatively coincide with the last closed magnetic surfaces (dotted lines in Figure 1(b)).

The above alignment suggests the tendency of the system to produce force-free states, where $\mathbf{J} = (4\pi)^{-1}\nabla \times \mathbf{B} = \lambda\mathbf{B}$, as in large-scale fluid models (Taylor 1974). In this regard, we show

the force-free parameter $\lambda(r) = \mathbf{J} \cdot \mathbf{B}/B^2$ in Figure 2(b). In contrast to the classical and global (constant-λ) force-free states, our local version of the minimization process (for each vortex) reveals a strong dependency of λ as a function of r, approaching zero at the vortex boundaries—namely, a non-linear force-free state with a uniform twist per unit length. This radial dependency of λ might suggest a more complex relaxation process (Montgomery et al. 1992; Servidio et al. 2010b).

### 4. A Kinetic Model

In what follows, we discuss how to interpret the long-lived structures in terms of kinetic plasma theory. We recall that the Vlasov equation for the αth species distribution function $f_\alpha(\mathbf{x}, \mathbf{v}, t)$ can be written as $\partial f_\alpha/\partial t = \{H, f_\alpha\}$, where H is the particle Hamiltonian and $\{\cdot, \cdot\}$ are the standard Poisson brackets. A stationary equilibrium is thus characterized by $\{H, f_\alpha\} = 0$, where $f_\alpha$ must be represented as a function of the integrals of motion. Inspired by the popular Harris approach (Harris 1962), we consider an exponential dependence on the invariants (energy and momenta), neglecting relativistic corrections (the bulk flows are nonrelativistic), with a simple drifting Maxwellian in a cylindrical coordinate system,

$$f_{\alpha}(\mathbf{x}, \mathbf{v}) = f_{\alpha 0} \exp \left[ - \frac{\mathcal{E}_{\alpha} - v_{\alpha}^{*} P_{z\alpha} - \Omega_{\alpha}^{*} P_{\phi\alpha}}{k_{B} T_{\alpha}} \right], \quad (1)$$

where $f_{\alpha 0} := N_{\alpha 0}(m_\alpha/2\pi k_B T_\alpha)^{3/2}$ is a normalization constant, $T_\alpha$ is the temperature, $\mathcal{E}_\alpha := m_\alpha(v_r^2 + v_\phi^2 + v_z^2)/2 + q_\alpha \psi(r)$ is the particle energy, and $P_{\phi\alpha} := r[m_\alpha v_\phi + q_\alpha A_\phi(r)/c]$ and $P_{z\alpha} := m_\alpha v_z + q_\alpha A_z(r)/c$ are the azimuthal and vertical momentum, respectively ($A_\alpha$ are the components of the averaged vector potential in the Lorenz gauge, $\psi$ is the electrostatic potential, and $v_\alpha$ are the components of the particle velocities). As usual in the literature, $k_B$ is the Boltzmann constant, $c$ is the speed of light, and $m_\alpha$ is the rest mass of each species. The undetermined (free) quantities appearing in Equation (1) stand for the out-of-plane characteristic linear velocity $v_\alpha^*$ and azimuthal velocity $\Omega_\alpha^*$. Once these quantities are specified, an exact kinetic equilibrium can be constructed.

By taking the moments of Equation (1), it is possible to compute the particle number densities $N_\alpha(r) = \int f_\alpha d^3\mathbf{v} = \exp[\gamma_\alpha(r)]$, where

$$\gamma_{\alpha}(r) = \ln(N_{0\alpha}) + \frac{q_{\alpha}}{c k_{B} T_{\alpha}} [v_{\alpha}^{*} A_{z}(r) + r \Omega_{\alpha}^{*} A_{\phi}(r) - c \psi(r)] + \frac{m_{\alpha}}{2 k_{B} T_{\alpha}} (r^{2} \Omega_{\alpha}^{*2} + v_{\alpha}^{*2}). \quad (2)$$

Similarly, by computing the azimuthally averaged bulk velocity $U(r) = N^{-1}(r)\int vfd^3v$, it is simple to demonstrate that it is related to the free parameters, being $U_{\alpha z}(r) = v_{\alpha z}^*$ and $U_{\alpha \phi}(r) = \Omega_{\alpha}^* r$. The problem can be further simplified with an assumption similar to the one made within the Harris sheet pinch, i.e., by imposing a (negligible) net constant charge
<figure><img src="image_7.png" /><figcaption>Figure 3. (a) Total magnetic field and current density components for the central vortex in Figure 1 (open symbols) and the KVR model (solid lines). Data have been averaged over time in the vortex frame.</figcaption></figure>

density $\kappa := \exp[\gamma_e(r)] - \exp[\gamma_p(r)]$ that yields

$$ 
\begin{aligned}
v_p^* &= -\frac{T_p}{T_e}v_e^*, \quad \Omega_p^* = -\frac{T_p}{T_e}\Omega_e^*, \\
\kappa &= \frac{m_p\Omega_p^{*2}T_e - m_e\Omega_e^{*2}T_p}{2\pi e^2(T_p + T_e)}.
\end{aligned} 
\quad (3) 
$$

By using Equation (1) into Maxwell's equations via the densities and the current expressions, assuming stationarity $(\partial/\partial t=0)$, and imposing the simplifications of Equation (3), it is possible to obtain a set of ordinary differential equations (ODEs) constituting our kinetic vortex reconstruction (KVR) model:

$$ \frac{dA_z}{dr} = -B_{\phi}(r), \quad (4) $$

$$ \frac{dA_{\phi}}{dr} = -\frac{A_{\phi}(r)}{r} + B'_{z}(r), \quad (5) $$

$$ \frac{dB'_{z}}{dr} = \frac{4\pi e}{c}(\Omega_{e}^{*} - \Omega_{p}^{*})r e^{\gamma_{e}(r)} + \Omega_{p}^{*}\kappa r, \quad (6) $$

$$ \frac{dB_{\phi}}{dr} = -\frac{B_{\phi}(r)}{r} + \frac{4\pi e}{c}(v_{p}^{*} - v_{e}^{*})e^{\gamma_{e}(r)} - v_{p}^{*}\kappa. \quad (7) $$

Here, $B'_z$ denotes the fluctuations in the out-of-plane field, i.e., $B_z = B'_z + b_{0z}$. Equations (4)–(7) can be integrated numerically after specifying, for each vortex, the internal temperatures of the eye $T_\alpha$ and the free parameters $(v_e^*, \Omega_e^*)$. The latter are obtained through a data-driven Monte Carlo method. A discussion of the reconstruction process can be found in Appendices B and C. It should be emphasized that a charge separation, albeit small (less than 10%), is present, especially in the vortex eye (not shown here).

Figure 3 offers a direct comparison between the KVR model (lines) and the actual numerical data (symbols). The model captures the behavior of all the fields and does so particularly well in the inner regions of the vortex ($r \lesssim r_0$) and for the magnetic field dependence. Note the presence of an azimuthal

<figure><img src="image_8.png" /><figcaption>Figure 4. Merging history of two long-lived structures as a function of Alfvénic crossing time (from top to bottom). The plot shows the plasma density (left column), the (force-free) cosine angle (middle column), and the magnetic helicity density (right column). Circumferences (in black) represent the eye of each vortex.</figcaption></figure>

current $J_\phi$ vanishing at the center and asymptotically and a nonzero vertical current $J_z$ at the axis.

Even outside the eye, the magnetic field components are qualitatively consistent with the model, even though the data are subject to the cyclonic arms, observed in Figures 1 and 2. Very similar behaviors are also shown by the other vortices tracked during the evolution: $B_z$ manifests a peak in the eye and diminishes for large $r$, whereas $B_\phi$ peaks at the vortex wall.

To understand the merging dynamics of existing vortices, we examine the simulation at different times, identifying and monitoring the most stable and long-lasting structures. In particular, we report in Figure 4 the “merging history” of two individual, isolated vortices (or magnetic islands) that wander across the background until they eventually merge with a companion. The figure shows, as a function of time (from top to bottom), the plasma number density, the cosine angle, and the magnetic helicity density (see later). If the structures are both very energetic, the encounter is quite explosive, resulting in a net current layer between them, where magnetic reconnection occurs (Servidio et al. 2009; Ripperda et al. 2019) and nonthermal particles are produced (Comisso & Sironi 2019). The $\cos(j, b)$ is quite strong inside the structures, and it changes sign in between them. The example reported in Figure 4 also helps us to appreciate that the spiral arms described above are actually the heritage of merger events. The process is related to the conservation of magnetic helicity, in a picture similar to Alexakis et al. (2006), where the magnetic helicity inverse cascade can occur at all scales. In our case, this results in metastable vortices ranging from large injection scales (a few $d_p$) to electron scales within the subinertial range.

We now reconcile the KVR model in Figure 3 with the fluidlike states in Figure 2. By taking the moments of the Vlasov equation, it is possible to define a hierarchy of the different contributions intervening in determining the equilibrium in the islands and conclude that the magnetic forces dominate the dynamics of the long-lived equilibria (see Appendix D for details). As a result, neglecting secondary effects such as the charge separation, one can combine
Equations (6) and (7) as

$$ \nabla \times \mathbf{B} = \frac{4\pi}{c} \mathbf{J}_0^* e^{\gamma_e(r)}, \qquad (8) $$

where $\mathbf{J}_0^* := q_e[0, (\Omega_p^* - \Omega_e^*)r, (v_p^* - v_e^*)]$ is a current-to-number density and $\gamma_e(r)$ is defined in Equation (2). When $r \to 0$, all functions, including $\gamma_e(r)$ and $\psi(r)$, are regular, smooth, and continuous (see Appendix B). Furthermore, for all vortices, we observe that $B_z$ peaks near the origin at a nonzero value, and the azimuthal component $B_\phi \sim \alpha r$, where $\alpha$ is a constant. Phenomenologically, to suppress the Lorentz force term in the momentum equation, the magnetic field $\mathbf{B}$ tends to align with $\mathbf{J} \approx \mathbf{J}_0^* e^{\gamma_e(r)}$, so Equation (8) can be approximated as

$$ \nabla \times \mathbf{B} = f(r) \mathbf{J}_0^* \sim \lambda(r) \mathbf{B}(r), \qquad (9) $$

where all the rescaling constraints are included in $\lambda(r)$. A very robust solution to Equation (9) is known as the Gold-Hoyle (GH) vortex (Gold & Hoyle 1960)—a flux tube describing a force-free, twisted field. This equilibrium was discovered for force-free coronal structures and has potential uses in astrophysical contexts (Fushiki & Sakai 1995). It is given by Allanson et al. (2016),

$$ A(\tilde{r}) = \frac{B_0}{2\xi} \left(0, \frac{1}{\tilde{r}} \ln(1 + \tilde{r}^2), \mp \ln(1 + \tilde{r}^2)\right), \qquad (10) $$

$$ \mathbf{B}(\tilde{r}) = B_0 \left(0, \pm \frac{\tilde{r}}{1 + \tilde{r}^2}, \frac{1}{1 + \tilde{r}^2}\right), \qquad (11) $$

$$ \mathbf{J}(\tilde{r}) = 2 \frac{\xi B_0}{\mu_0} \left(0, \frac{\tilde{r}}{(1 + \tilde{r}^2)^2}, \pm \frac{1}{(1 + \tilde{r}^2)^2}\right), \qquad (12) $$

where $B_0$ is a typical field strength, $\tilde{r} = \xi r$ is a dimensionless radial coordinate, and $\xi$ is a characteristic gradient. The solution is valid for “strong” vortices (fluctuations of the order of the mean field) and generalized for clockwise (upper sign) or counterclockwise (lower sign) rotation.

To determine whether the above model is a universal property of the turbulent cascade, we detect and analyze all the ten long-lived structures appearing during the simulation, collected in Figure 5. Here, the field components are normalized to $B_0$ and rescaled via $\xi$, which is determined through a simple fitting procedure; also reported are the corresponding GH magnetic components from Equation (11) (to avoid overcrowding, we represent $-|B_\phi|$). The GH solution provides a qualitatively accurate description, implying that the KVR (and its GH approximation) may characterize energetic structures in plasma turbulence. This kind of solution is achieved through a local relaxation process where the magnetic helicity is finite and plays a crucial role (Woltjer 1958; Taylor 1974; Matthaeus et al. 1982; Alexakis et al. 2006). Such an MHD invariant measures the twisting of the field lines and, in 2.5D, is defined as the volume average $H_m = V^{-1} \int_V h_m d^3x$, where $h_m = \mathbf{a}_\perp \cdot \mathbf{b}_\perp$, with $\mathbf{a}_\perp = (a_x, a_y)$ and $\mathbf{b}_\perp = (b_x, b_y)$ being the in-plane components. Contrary to the classical Bessel solutions of the linear force-free state, the numerical results are more consistent with a constant twist per unit length (along the vortex axis). As shown in Figure 4, the structures retain a finite amount of magnetic helicity before and after the merging. From Equations (10) and (11), performing a volume average over the flux tube, one gets $H_m \sim 10^{-1} B_0^2 r_0$ (e.g., $B_0 \approx 43$ and $r_0 \approx 170 d_e$ for vortex 1 in Figure 1(b)), which is in accordance

<figure><img src="image_9.png" /><figcaption>Figure 5. Radial behavior of the magnetic field components of all the long-lived vortices appearing during the simulation (points). Components have been rescaled to the field at the center $B_0$ and using their typical gradient $\xi$. The GH solution (shaded) describes qualitatively well the profiles near the eye.</figcaption></figure>

with the observed values in Figure 4 ($H_m \sim 3 \times 10^4$). It is particularly interesting to note that all the structures described here in the context of kinetic theory closely resemble those observed in the magnetosheath (Alexandrova et al. 2006) and the solar wind, with sizes spanning from MHD to sub-ion scales (Vinogradov et al. 2023).

## 5. Conclusions

Exploiting the results of direct numerical PIC simulations, we propose a description of plasma turbulence envisioned as a mosaic of equilibrium-like patterns. In such a scenario, coherent structures emerge from the turbulent background and occasionally encounter and merge with other similar metastable structures during their life. This self-similar process systematically produces newborn equilibria, obeying a kinetic stationary solution of the Vlasov equation. The magnetic vortices qualitatively obey a universal form that can be simply characterized using the GH equilibrium. These structures show a characteristic size typical of inertial range turbulence, thus making them macroscopically relevant: they might grow, by coalescence, to a significant fraction of astrophysical system size, with potential observable signatures (El Mellah et al. 2022; Vos et al. 2024). The present work focuses on a 2.5D model, considering the plane perpendicular to a mean magnetic field. While this 2D approximation differs qualitatively from the (more complex and expensive) 3D case, it may still provide insights into some relaxation processes characteristic of magnetized astrophysical plasmas. Indeed, these equilibria might also be relevant for full 3D anisotropic settings, in cases where an external field effectively reduces the dimensionality of turbulence (Shebalin et al. 1983; Chernoglazov et al. 2021; Khabarova et al. 2021; Ripperda et al. 2022). In such a general case, which will motivate future investigations, the KVR model can acquire a weak dependency along the magnetic field coordinate, say $z$, as typical of solar flux ropes. We plan to investigate the dynamics of magnetic helicity in these vortices on macroscopic scales (in a fluidlike regime) by using larger domains. Additionally, in future work, we will describe the
formation process of each individual vortex and how the GH solution is achieved over time, while also exploring the plasma parameter space.

Since the long-lived structures might potentially grow, our results could be significantly relevant for the comprehension of astrophysical plasmas, especially in scenarios where transient and flare emissions are associated with the formation of plasmons during the accretion process, as observed, for instance, in Sgr A*, as well as in the observation of flux tubes in the solar wind and corona (varying the magnetization $\sigma$ and the plasma $\beta$ parameters).

## Acknowledgments

The authors thank Luciano Rezzolla and William H. Matthaeus for useful discussions. C.M. acknowledges the support from the ERC Advanced Grant “JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales” (grant No. 884631). B.C. acknowledges the support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 863412). S.S. acknowledges “Progetto STAR 2-PIR01 00008” (Italian Ministry of University and Research). A.C.O. gratefully acknowledges the “Ciencia Básica y de Frontera 2023–2024” program of the “Consejo Nacional de Humanidades, Ciencias y Tecnologías” (CONAHCYT, Mexico), projects CBF2023-2024-1102 and 257435. Computational resources were provided by CINECA through the ISCRA Class B project “KITCOM—HP10BB7U73.” Finally, the authors would like to thank the anonymous reviewer for the very useful comments and suggestions.

## Appendix A
The Numerical Method

Simulations are performed by means of the PIC code Zeltron, which solves the following system of equations (Lorentz–Newton + Maxwell’s equations) for particles and field evolution:

$$v_i = \frac{dr_i}{dt} = \frac{cu_i}{W_i}, \qquad (A1)$$

$$\frac{du_i}{dt} = \frac{q_\alpha}{m_\alpha c} \left[ \boldsymbol{\epsilon} + \frac{v_i \times \boldsymbol{b}}{c} \right], \qquad (A2)$$

$$\nabla \cdot \boldsymbol{\epsilon} = 4\pi\rho_c, \qquad (A3)$$

$$\nabla \cdot \boldsymbol{b} = 0, \qquad (A4)$$

$$\frac{\partial \boldsymbol{b}}{\partial t} = -c \nabla \times \boldsymbol{\epsilon}, \qquad (A5)$$

$$\frac{\partial \boldsymbol{\epsilon}}{\partial t} = c \nabla \times \boldsymbol{b} - \mu_0 \boldsymbol{j}. \qquad (A6)$$

The first two equations describe the motion of particles in the Lagrangian specification and constitute the characteristic curves along which the Vlasov equation can be solved. Here, $\mathbf{r}_i$ and $\mathbf{v}_i$ are the position and the proper 3-velocity of the $i$th macroparticle (Cerutti et al. 2013), $\mathbf{u}_i = \mathbf{p}_i (m_\alpha c)^{-1}$ represents its normalized momentum as measured by an inertial observer at rest ($\mathbf{p}_i$ being the particle’s momentum in that frame), and $W_i = (1 - v_i^2/c^2)^{-1/2}$ is the associated Lorentz factor. The subscript $\alpha$ stands for either protons ($p$) or electrons ($e$).

Regarding Maxwell’s equations, all quantities are well described at the beginning of Section 2 of this Letter.

We adopt the geometrized unit system, where the speed of light $c$, the gravitational constant $G$, the elementary charge $q_e$, the electron mass $m_e$, the Boltzmann constant $k_B$, and the reduced Planck constant $\hbar$ are set to unity, whereas the vacuum permittivity $\epsilon_0 = 1/4\pi$ and the vacuum permeability $\mu_0 = 4\pi$.

The code utilizes the Yee algorithm (Yee 1966) to solve the time-dependent Maxwell's equations, in which the different components of the fields are staggered in both space (on the grid) and time. Zeltron has a second-order error in space and time, ensuring that the magnetic constraint $\nabla \cdot \boldsymbol{b} = 0$ is satisfied at any time step of the simulation. All the quantities are expressed in terms of electron skin depths, i.e., $d_e := c/\omega_{pe} = c\sqrt{m_e/4\pi n_0 q_e^2} = 1$, being the number density at equilibrium $n_0 = n_e = n_p = (4\pi)^{-1}$.

The setup settings of the simulation are described in the following. We use periodic boundary conditions and a spatial resolution such that the mesh is $\Delta x = \Delta v = d_e/3$, employing 10 particles per cell, namely, five ions and five electrons, with the full mass ratio. Time resolution holds $\Delta t = 0.45 \Delta x = 0.15 d_e$. The magnetization (which determines the available magnetic energy per particle) $\sigma := b_{0z}^2/4\pi w_0$, where the out-of-plane mean magnetic field magnitude $b_{0z} \approx 43$ and the enthalpy density $w_0 := (n_p m_p + n_e m_e)c^2 + \Gamma_e \theta_e + \Gamma_p \varrho_p$, with the internal energy density $\varrho_\alpha \approx n_\alpha k_B T_\alpha (\Gamma_\alpha - 1)^{-1}$ and the adiabatic index $\Gamma_\alpha = 4/3$. Since we are in a weakly relativistic regime, the magnetization $\sigma$ is approximately equal to $\sigma_p^{\text{cold}} := b_{0z}^2/4\pi n_p m_p c^2$. The choice $\sigma \approx \sigma_p^{\text{cold}} = 1$ indicates an equilibrium between magnetic and kinetic forces, which is typical in astrophysical scenarios such as black hole winds near accretion disks. Additionally, the plasma beta is $\beta_p = \beta_e := 8\pi n_0 k_B T_\alpha / b_{0z}^2 = 3 \times 10^{-3}$, where the temperature $T_\alpha := \theta_\alpha m_\alpha c^2 / k_B$ is chosen to ensure initial thermal equilibrium between species. We note that the choice of both $\sigma$ and $\beta$ is crucial for determining the regime of plasma (Ball et al. 2018). The adimensional temperatures for each species are equal to $\theta_p = 1.5 \times 10^{-3}$ and $\theta_e \approx 2.75$, respectively. Such a configuration allows us to largely resolve the Debye plasma length $\lambda_D := \sqrt{(m_e \theta_e + m_p \theta_p)c^2/4\pi n_0 e^2} \approx 2.35 d_e$. Besides, when the turbulence is fully developed, the velocity distribution of electrons is highly nonthermal, and their Larmor radius becomes significantly larger due to huge accelerations, effectively increasing our resolution. We further define the Alfvénic crossing time as $t_A := L_0/v_A \approx 7748$; hence, the Alfvén velocity reads $v_A := c\sqrt{\sigma/1 + \sigma} \approx 0.71$.

We impose large-scale, random initial conditions for the magnetic field power spectrum to achieve a strong turbulent state. To avoid excessive compressive activity, no out-of-plane magnetic field fluctuations are prescribed at the beginning, and no bulk flows or density perturbations are initiated.

The analysis for our fiducial simulation is carried out at a time $t \approx 3.23 t_A$ when the peak of nonlinear activity has long been reached. Indeed, as shown in Figure 6, the rms of the vertical component of the current density $\mathbf{j}_z = (\nabla \times \boldsymbol{b}/4\pi) \cdot \hat{\mathbf{z}}$ experiences an absolute maximum much earlier, around $t \approx 0.62 t_A$.
<figure><img src="image_10.png" /><figcaption>Figure 6. The rms of the vertical component of the current density (black solid line). Dotted red and dashed blue lines represent the time of most intense nonlinear activity and the analysis time, respectively.</figcaption></figure>

## Appendix B
### The Kinetic Equilibrium Principle

As discussed in this Letter, in line with the Harris equilibrium, we consider a particle distribution function with an exponential dependence on the three invariants, as represented in Equation (1). Following this ansatz, the probability density of finding a particle at a given point in the phase space can be explicitly written as

$$f_{\alpha}(r, v_r, v_{\phi}, v_z) = N_{0\alpha} \left( \frac{m_{\alpha}}{2\pi k_{\text{B}} T_{\alpha}} \right)^{3/2} \exp \left[ -\frac{m_{\alpha}}{2k_{\text{B}}T_{\alpha}} (v_r^2 + v_{\phi}^2 + v_z^2) + \frac{m_{\alpha}v_{\alpha}^*}{k_{\text{B}}T_{\alpha}} v_z + \frac{m_{\alpha}\Omega_{\alpha}^*}{k_{\text{B}}T_{\alpha}} r v_{\phi} + \frac{q_{\alpha}v_{\alpha}^*}{c k_{\text{B}}T_{\alpha}} A_z(r) + \frac{q_{\alpha}\Omega_{\alpha}^*}{c k_{\text{B}}T_{\alpha}} r A_{\phi}(r) - \frac{q_{\alpha}}{k_{\text{B}}T_{\alpha}} \psi(r) \right], \quad (B1)$$

which must satisfy Maxwell's equations:

$$\begin{align*}
\nabla \cdot \mathbf{E} &= 4\pi e \left( \int_{\Omega_v} f_p(r, \mathbf{v}) d^3\mathbf{v} - \int_{\Omega_v} f_e(r, \mathbf{v}) d^3\mathbf{v} \right), \\
\nabla \cdot \mathbf{B} &= 0, \\
\nabla \times \mathbf{E} &= -\frac{1}{c} \frac{\partial \mathbf{B}}{\partial t} = \mathbf{0}, \\
\nabla \times \mathbf{B} &= \frac{4\pi e}{c} \left( \int_{\Omega_v} \mathbf{v} f_p(r, \mathbf{v}) d^3\mathbf{v} - \int_{\Omega_v} \mathbf{v} f_e(r, \mathbf{v}) d^3\mathbf{v} \right). \tag{B2}
\end{align*}$$

We point out that we are assuming stationarity, i.e., $\partial/\partial t = 0$. Given that the fields depend only on the radial coordinate, thus $\partial/\partial\phi = \partial/\partial z = 0$, Equations (B2) reduce to three ODEs for the potentials (one for the electric potential and two for the

azimuthal and vertical components of the vector potential):

$$\begin{align*}
\frac{1}{r} \frac{d}{dr} \left( r \frac{d\psi}{dr} \right) &= -4\pi e \left( \int_{\Omega_v} f_p d^3\mathbf{v} - \int_{\Omega_v} f_e d^3\mathbf{v} \right), \\
\frac{d}{dr} \left( \frac{1}{r} \frac{d(rA_{\phi})}{dr} \right) &= -\frac{4\pi e}{c} \left( \int_{\Omega_v} v_{\phi} f_p d^3\mathbf{v} - \int_{\Omega_v} v_{\phi} f_e d^3\mathbf{v} \right), \\
\frac{1}{r} \frac{d}{dr} \left( r \frac{dA_z}{dr} \right) &= -\frac{4\pi e}{c} \left( \int_{\Omega_v} v_z f_p d^3\mathbf{v} - \int_{\Omega_v} v_z f_e d^3\mathbf{v} \right). \tag{B3}
\end{align*}$$

By taking the moments of Equation (1), it is possible to compute the particle density $N_{\alpha}(r) = \exp[\gamma_{\alpha}(r)]$, where $\gamma_{\alpha}(r)$ is given by Equation (2), and the bulk velocities $U_{\alpha z}(r) = v_{\alpha z}^*$ and $U_{\alpha \phi}(r) = \Omega_{\alpha}^* r$. The physical meaning of the free parameters $v_e^*$ and $\Omega_e^*$ is specified in Section 4 of the leading part of the current work, where the kinetic model is illustrated. Having these quantities, the current densities $J_{\phi \alpha}(r) = q_{\alpha} N_{\alpha}(r) U_{\phi \alpha}(r)$ and $J_{z \alpha}(r) = q_{\alpha} N_{\alpha}(r) U_{z \alpha}(r)$ are derived straightforwardly.

We observe a small and approximately constant net charge separation, which we define as

$$\kappa = N_e(r) - N_i(r) = \exp[\gamma_e(r)] - \exp[\gamma_p(r)], \quad (B4)$$

where $\kappa$ is a constant; hence, the calculation can be further simplified by assuming $\nabla \cdot \mathbf{E} = -\nabla^2 \psi = -4\pi e \kappa$. This assumption leads directly to a straightforward choice for the electric potential of the form $\psi(r) = \psi(0) + (\pi e \kappa) r^2$. Due to gauge invariance, we set all potentials to zero at the vortex axis, i.e., $\psi(0) = 0$, $A_z(0) = 0$, and $A_{\phi}(0) = 0$. Introducing the quantity $\eta(r) = \gamma_e(r) - \gamma_p(r)$, Equation (B4) becomes

$$\kappa = \exp[\gamma_e(r)] - \exp[\gamma_e(r) - \eta(r)], \quad (B5)$$

which immediately yields $\eta(r) = -\ln\{1 - \kappa \exp[-\gamma_e(r)]\}$. The relationship between the exponents $\gamma_e(r)$ and $\gamma_p(r)$ leads to

$$\begin{align}
C_p + \frac{q_p v_p^*}{c k_B T_p} A_z(r) + \frac{q_p \Omega_p^*}{c k_B T_p} r A_{\phi}(r) - \frac{q_p}{k_B T_p} \psi(r) & \nonumber \\
+ \frac{m_p}{2 k_B T_p} \Omega_p^{*2} r^2 + \frac{m_p}{2 k_B T_p} v_p^{*2} & = C_e + \frac{q_e v_e^*}{c k_B T_e} A_z(r) \nonumber \\
+ \frac{q_e \Omega_e^*}{c k_B T_e} r A_{\phi}(r) - \frac{q_e}{k_B T_e} \psi(r) + \frac{m_e}{2 k_B T_e} \Omega_e^{*2} r^2 & \nonumber \\
+ \frac{m_e}{2 k_B T_e} v_e^{*2} - \eta(r) & \tag{B6}
\end{align}$$

where $\eta(r)$ is a small function due to the minimal net charge separation (consistent with quasi-neutrality) and $C_{\alpha} := \ln(N_{0\alpha})$. By employing the correspondence between the exponents at the vortex center, where $\gamma_p(0) = \gamma_e(0) - \eta(0)$, we derive clear expressions for $C_p$ and $C_e$, which are then introduced into Equation (B6). Noticing that $\eta(r) - \eta(0)$ is a negligible term, we perform some algebra to obtain the following polynomial relation:

$$\zeta_0 \psi(r) + \zeta_1 A_z(r) + \zeta_2 A_{\phi}(r) r + \zeta_3 r^2 = 0. \quad (B7)$$

A possible solution can be found by setting $\zeta_1 = 0$ and $\zeta_2 = 0$, from which the conditions shown in Equation (3) follow. The latter recalls the Harris equilibrium in the sheet
pinch for a nonhomogeneous magnetic field in cylindrical geometry, with an additional condition for the in-plane azimuthal flow. Under this choice, a simple relation for the electric potential can be made explicit:

$$ \psi(r) = -\frac{\zeta_3}{\zeta_0} r^2 = (\pi e \kappa) r^2 = \frac{m_p \Omega_p^{*2} T_e - m_e \Omega_e^{*2} T_p}{2e(T_p + T_e)} r^2. \quad (B8) $$

After that, we invoke Equation (B1) in the last two differential equations of the system (Equation (B3)) to finally arrive at the set of governing Equations (4)–(7) of the main Letter, which represents the core of the KVR model.

## Appendix C
Data-driven Optimization Model

The KVR system of equations can be solved through direct numerical integration using a second-order Runge–Kutta scheme once the typical bulk electron velocity $v_e^*$ and angular velocity $\Omega_e^*$ are given, and boundary conditions at each vortex center are applied. For each coherent structure, we search for this pair of free parameters by implementing a data-driven Monte Carlo numerical technique capable of selecting the best fit for each vortex, starting from the results of the integration on a parameter space $S$ given by $\{v_e^*, \Omega_e^*\}$.

Inspecting the most representative vortex in Figure 1(b) at a time beyond the maximum of nonlinear activity ($t \approx 3.23t_A$), we impose $B_z'(0) \approx 40.5$, $B_\phi(r)/r|_{r=0} = 0$, and $A_\phi(r)/r|_{r=0} = 0$ as boundary conditions. As for the species temperatures ($T_e$, $T_p$), we choose the averaged parallel temperature profiles in the vortex eye, i.e., $T_\alpha = T_{\parallel\alpha}$, with $T_{\parallel e} = 51.8$ and $T_{\parallel p} = 51.7$.

Each long-lived structure is uniquely determined by the choice of the star parameters $v_e^*$ and $\Omega_e^*$. Such a choice has to ensure that the distribution function, once integrated over the phase space $\Omega$, returns the different fields and momenta in the best possible agreement with the data. Hence, a Monte Carlo method is performed over the parameter space $S$, calculating the discrepancies between the KVR outputs and the observed profiles as below.

We span the parameter space by discretizing the domain $\{v_e^*, \Omega_e^*\} \in [-0.2, 0.2], [-0.0015, 0.0015]$, performing a campaign of reconstructions to minimize the error

$$ \varepsilon_g(v_e^*, \Omega_e^*) = \frac{\int_0^{r_0} [g_{\text{KVR}}(r) - g_{\text{sim}}(r)]^2 dr}{\int_0^{r_0} g_{\text{sim}}(r)^2 dr}, \quad (C1) $$

where $g_{\text{sim}}$ is a generic field from the simulation and $g_{\text{KVR}}$ is the reconstructed field from Equations (4)–(7). We minimize for $g = [A_\phi, A_z, B_\phi, B_z', J_\phi, J_z, N_e, N_p]$, attaining for each of them the associated error $\varepsilon_g$. To achieve the best possible profiles, we estimate a single average over all the fields, namely, $\varepsilon_{\text{avg}}(v_e^*, \Omega_e^*) = \bar{\varepsilon}_g(v_e^*, \Omega_e^*)$.

In Figure 7, we show a zoom of the average error $\varepsilon_{\text{avg}}$, with its minimum located in the fourth quadrant of the parameter space (positive $v_e^*$ and negative $\Omega_e^*$). The minimum is placed at $v_e^* = 0.09$ and $\Omega_e^* = -5.30 \times 10^{-4}$, as indicated by a cyan circle.

The resulting parameters are closely related to the flux tube current density: for the KVR model, we get $\mathbf{J}(r) := q_e[0, (N_p(r)\Omega_p^*-N_e(r)\Omega_e^*)r, (N_p(r)v_p^*-N_e(r)v_e^*)]$, from which $\langle J_\phi^* \rangle_r \approx 1.17 \times 10^{-2}$ and $\langle J_z^* \rangle_r \approx -2.73 \times 10^{-2}$; for the observed structures in Figure 1(b), we measure

<figure><img src="image_11.png" /><figcaption>Figure 7. Average error estimated by considering all the vortex profiles. The best choice for the free parameters is highlighted by a cyan circle.</figcaption></figure>

$\langle J_\phi \rangle_r \approx 1.27 \times 10^{-2}$ and $\langle J_z \rangle_r \approx -4.69 \times 10^{-2}$. Such an estimate indicates a very good agreement between the model and the data.

## Appendix D
Momentum Equation and Hierarchy of Forces

Computing the moments of the Vlasov–Maxwell system, one finds the momentum equation, which, in a stationary case, reads

$$ \nabla \cdot (\rho_m \mathbf{U} \otimes \mathbf{U}) + \nabla \cdot \mathbf{P}_{\text{tot}} - \rho_c \mathbf{E} - \mathbf{J} \times \mathbf{B} = \mathbf{0}, \quad (D1) $$

where $\rho_m$ represents the total mass density, $\rho_c$ is the total charge density, and $\mathbf{P}_{\text{tot}}$ is the total pressure tensor. Employing the continuity equation in the above expression, recalling the divergence of a symmetric tensor of rank two in cylindrical coordinates, and imposing $U_r = 0$, $\partial/\partial\phi = \partial/\partial z = 0$, one can decompose Equation (D1) along the three cylindrical axes:

$$ \hat{r}: -\frac{\rho_m U_\phi U_\phi}{r} + \left(\frac{1}{r}\partial_r(r P_{rr}) - \frac{P_{\phi\phi}}{r}\right) - \rho_c E_r - (J_\phi B_z - J_z B_\phi) = 0, \quad (D2) $$

$$ \hat{\phi}: \left( \partial_r P_{r\phi} + \frac{2P_{r\phi}}{r} \right) - \rho_c E_{\phi} - (J_z B_r - J_r B_z) = 0, \quad (D3) $$

$$ \hat{z}: \frac{1}{r} \partial_r (r P_{rz}) - \rho_c E_z - (J_r B_\phi - J_\phi B_r) = 0. \quad (D4) $$

Since the radial components of currents and magnetic fields are null, the out-of-diagonal terms of the pressure are negligible, and the main electric field is radial, it is worth noting that the only nontrivial equation is the radial one.

Establishing a hierarchical order for all the terms appearing in the radial momentum equation proves helpful in extracting the critical constituents of such metastable states. We compute all the quantities of Equation (D2) and compare them in Figure 8. Computing radial averages for each of them, it is evident that magnetic forces dominate the dynamics; in fact, a clear ordering can be inferred: $|J_i B_j|_{i,j \neq r} > |\rho_c E_r| \gg |(\nabla \cdot \mathbf{P}_{\text{tot}}) \cdot \hat{r}| > |[\nabla \cdot (\rho_m \mathbf{U} \otimes \mathbf{U})] \cdot \hat{r}|$.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Components of the momentum equation for the fluid approximation in the vortex eye.</figcaption></figure>

## ORCID iDs

Mario Imbrogno https://orcid.org/0000-0002-7459-5735

Claudio Meringolo https://orcid.org/0000-0001-8694-3058

Sergio Servidio https://orcid.org/0000-0001-8184-2151

Alejandro Cruz-Osorio https://orcid.org/0000-0002-3945-6342

Benoît Cerutti https://orcid.org/0000-0001-6295-596X

Francesco Pegoraro https://orcid.org/0000-0002-7216-5491

## References

* Alexakis, A., Mininni, P. D., & Pouquet, A. 2006, ApJ, 640, 335

* Alexandrova, O. 2008, NPGeo, 15, 95

* Alexandrova, O., Mangeney, A., Maksimovic, M., et al. 2006, JGRA, 111, 12208

* Alexandrova, O., Saur, J., Lacombe, C., et al. 2009, PhRvL, 103, 165003

* Allanson, O., Wilson, F., & Neukirch, T. 2016, PhPl, 23, 092106

* Baiotti, L., & Rezzolla, L. 2017, RPPh, 80, 096901

* Bale, S. D., Kellogg, P. J., Mozer, F. S., Horbury, T. S., & Reme, H. 2005, PhRvL, 94, 215002

* Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80

* Borovsky, J. E. 2008, JGRA, 113, 08110

* Carnevale, G. F., McWilliams, J. C., Pomeau, Y., Weiss, J. B., & Young, W. R. 1991, PhRvL, 66, 2735

* Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147

* Chavanis, P., & Sommeria, J. 1998, JFM, 356, 259

* Chernoglazov, A., Ripperda, B., & Philippov, A. 2021, ApJL, 923, L13

* Comisso, L., & Sironi, L. 2019, ApJ, 886, 122

* Drake, J., Opher, M., Swisdak, M., & Chamoun, J. 2010, ApJ, 709, 963

* El Mellah, I., Cerutti, B., Crinquand, B., & Parfrey, K. 2022, A&A, 663, A169

* Fushiki, T., & Sakai, J.-I. 1995, SoPh, 161, 317

* Gold, T., & Hoyle, F. 1960, MNRAS, 120, 89

* Goldstein, M. L., Roberts, D. A., & Matthaeus, W. 1995, ARA&A, 33, 283

* Haller, G. 2005, JFM, 525, 1

* Harris, E. G. 1962, NCim, 23, 115

* Karimabadi, H., Roytershteyn, V., Wan, M., et al. 2013, PhPl, 20, 012303

* Khabarova, O., Malandraki, O., Malova, H., et al. 2021, SSRv, 217, 38

* Lukin, A., Vasko, I., Artemyev, A., & Yushkov, E. 2018, PhPl, 25, 012906

* Matthaeus, W. H., Goldstein, M. L., & Smith, C. 1982, PhRvL, 48, 1256

* Matthaeus, W. H., Pouquet, A., Mininni, P. D., Dmitruk, P., & Breech, B. 2008, PhRvL, 100, 085003

* Matthaeus, W. H., Wan, M., Servidio, S., et al. 2015, RSPTA, 373, 20140154

* McWilliams, J. C. 1984, JFM, 146, 21

* Meringolo, C., Cruz-Osorio, A., Rezzolla, L., & Servidio, S. 2023, ApJ, 944, 122

* Montgomery, D., Matthaeus, W. H., Stribling, W. T., Martinez, D., & Oughton, S. 1992, PhFlA, 4, 3

* Nathanail, A., Mpisketzis, V., Porth, O., Fromm, C. M., & Rezzolla, L. 2022, MNRAS, 513, 4267

* Pecora, F., Greco, A., Hu, Q., et al. 2019, ApJL, 881, L11

* Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325

* Pokhotelov, O., & Petviashvili, V. 1992, Solitary Waves in Plasmas and in the Atmosphere (London: Gordon and Breach),

* Powell, M. D., & Houston, S. H. 1998, MWRv, 126, 1259

* Ripperda, B., Bacchini, F., & Philippov, A. A. 2020, ApJ, 900, 100

* Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJL, 924, L32

* Ripperda, B., Porth, O., Sironi, L., & Keppens, R. 2019, MNRAS, 485, 299

* Sahraoui, F., Goldstein, M. L., Robert, P., & Khotyaintsev, Y. V. 2009, PhRvL, 102, 231102

* Servidio, S., Matthaeus, W. H., & Dmitruk, P. 2008, PhRvL, 100, 095005

* Servidio, S., Matthaeus, W. H., Shay, M. A., Cassak, P. A., & Dmitruk, P. 2009, PhRvL, 102, 115003

* Servidio, S., Valentini, F., Califano, F., & Veltri, P. 2012, PhRvL, 108, 045001

* Servidio, S., Wan, M., Matthaeus, W. H., & Carbone, V. 2010a, PhFl, 17, 032315

* Servidio, S., Wan, M., Matthaeus, W. H., & Carbone, V. 2010b, PhFl, 22, 125107

* Shebalin, J. V., Matthaeus, W. H., & Montgomery, D. 1983, JPIPh, 29, 525

* Taylor, J. B. 1974, PhRvL, 33, 1139

* Vinogradov, A., Alexandrova, O., Démoulin, P., et al. 2023, arXiv:2307.10478

* Vos, J., Olivares, H., Cerutti, B., & Moscibrodzka, M. 2024, MNRAS, 531, 1554

* Woltjer, L. 1958, PNAS, 44, 833

* Yee, K. 1966, ITAP, 14, 302