Thermal equilibration in collisionless magnetospheric plasmas via entropy-mode turbulence

ABSTRACT

A planetary magnetosphere is a peculiar plasma environment where a high-temperature plasma is confined in a strongly inhomogeneous dipolar magnetic field generated by a planet. Turbulence driven by the magnetic curvature and density gradient (the entropy mode) is known to cause the inward pinch whereby particles are transported against the density gradient to achieve high confinement. However, a comprehensive understanding of how the global magnetospheric plasma confinement is determined is still missing. Here, we show that the entropy-mode turbulence equilibrates temperatures between species without collisions in a magnetospheric plasma. The classical stability analysis in terms of energetic consideration reveals the interchangeable roles of electrons and ions for destabilization due to resonance with drift waves depending on their temperature. One of the species has negative energy and grows its energy fed from the background density gradient. Turbulence driven by the microscopic (kinetic) instability tends to rearrange the internal energy between species, pre-dominantly via linear physics, leading to an equal temperature state. Our new finding adds an ingredient to energy transport processes and contributes to a consistent explanation of the global self-organization of magnetospheric plasmas.

instabilities, plasmas, turbulence

1 INTRODUCTION

A planetary magnetosphere is a peculiar plasma environment: a high-temperature plasma is confined in a strongly inhomogeneous dipolar magnetic field generated by a planet. Inspired by nature, the dipole confinement concept as a fusion device was proposed (Hasegawa 1987; Hasegawa, Chen & Mauel 1990) and later levitated ring dipole confinement devices have been built and operated for years (Boxer et al. 2010; Yoshida et al. 2010). It is ready to study this natural complex system in detail by in situ spacecraft observations (Burch et al. 2016; Nakamura et al. 2018; Benkhoff et al. 2021) and comprehensive diagnostics of experimental devices (Nishiura et al. 2017; Kenmochi et al. 2022; Saitoh et al. 2024) to gain insight into the interplay between the inhomogeneous field and plasma dynamics. Understanding the physics of turbulent self-organization processes occurring in magnetospheric plasmas opens a path to designing compact fusion machines favourable to plasmas which may be advantageous over the existing devices in which precise control of plasmas is necessary. In this paper, we report a newly discovered non-intuitive process that emerged from turbulence in the magnetosphere, that is the thermal equilibration without collisions, which may have a potential impact on the global self-organization of the magnetosphere. By examining the free energy driving turbulence, it is shown that the favourable state is a state where two species have equal temperatures by turbulent energy exchange between species.

The energy is a key to seeing how systems evolve. Thermodynamics tells us that closed systems tend to minimize the energy (with fixed entropy) toward thermodynamic equilibrium. While, if there are constraints, ordered structures appear (Hasegawa 1985). Zonal flows in two-dimensional turbulence (Rhines 1975; Hasegawa 1983; Hasegawa & Wakatani 1987) or the force-free state in magnetohydrodynamics (MHD) (Woltjer 1958; Taylor 1974) are the well-known examples. The free energy contained in non-equilibrium states drives turbulence to cause relaxation toward lower energy states. However, self-consistently determining the system-wide configuration of magnetospheric confinements is far beyond the scope of the present work since even possible kinetic equilibria are not at all known. Instead, we may ask, at given local states, in which direction the system evolves. Near the equilibrium, the instability leads to anomalous transport (temperature equilibration) if the free energy is contained in the background gradients of thermodynamic quantities (temperature difference between species). Such information will be utilized to construct a model to predict the density and temperature profiles of each species.

A prominent phenomenon induced by the inhomogeneity of magnetic configuration is the inward pinch whereby particles are transported against the density gradient, formally violating Fick’s law. The inward transport was primarily explained by the adiabatic theory of charged particle motion under fluctuating electric fields (Fälthammar 1965; Birmingham 1969). Later, the gyrokinetic theory applicable to slow collective dynamics was applied to analyse the linear stability and transport in the magnetospheric plasmas and revealed that the electrostatic entropy mode driven by the magnetic curvature and density gradient is responsible for turbulent transport (Simakov, Catto & Hastie 2001; Simakov, Hastie & Catto 2002). Gyrokinetic simulations of the entropy mode in the dipole configuration successfully demonstrated the existence of a particle pinch regime (Kobayashi, Rogers & Dorland 2009, 2010). Now, it is understood that the anomalous transport occurs to rearrange the profile by releasing free energy to achieve the marginally stable state against the MHD instability (Boxer et al. 2010; Kobayashi et al. 2010). An interesting explanation of this phenomenon in terms of the statistical mechanical point of view of Hamiltonian systems was proposed by Sato & Yoshida (2016). The up-hill diffusion results from entropy minimization in properly constructed constrained phase space incorporating the magnetic field structure.

Turbulence also heats plasmas (Howes et al. 2011; Barnes, Parra & Dorland 2012; Candy 2013; Barnes, Abiuso & Dorland 2018; Kawazura, Barnes & Schekochihin 2019; Kato et al. 2024). The turbulent heating can be a possible candidate to determine a thermal equilibrium state in collisionless magnetospheric environments, as well as astronomical or space objects, such as the accretion flows (Kawazura et al. 2019), solar winds (Howes et al. 2011), and fusion devices (Barnes et al.2012; Candy 2013; Kato et al. 2024). As will be shown later, the gyrokinetic formulation of the entropy mode in magnetospheric plasmas shows that electrons and ions behave in a similar manner independent of the mass ratio except for finite Larmor radius (FLR) effects. The role played by electrons and ions can be interchangeable. We address which of the species is responsible for driving the system in terms of free energy. It turns out that the negative energy due to wave–particle interactions carried by one species provides the free energy for destabilizing the system and the role played by each species changes according to its temperature. The free energy is of a kinetic origin. Therefore, the instability is supposed to relax the temperature difference by exchanging internal energy between species, as well as macroscopic rearrangement of profiles. The argument is essentially the same as that of Barnes et al. (2018) for the universal instability where the instability is caused by wave–particle interactions along the field line and therefore the thermal speed rather than the temperature itself matters.

In Section 2, we describe the electrostatic gyrokinetic model governing low-frequency dynamics in magnetospheres. We first perform a linear stability analysis to determine the instability mechanism in Section 3. Then, we show that the instability leads to thermal equilibration by non-linear energy exchange between species in Section 4. A conclusion is given in Section 5.

2 GYROKINETIC MODEL

We consider low-frequency electrostatic fluctuations in the equatorial plane of a magnetosphere shown in Fig. 1. The fluctuations are characterized by drift-type modes driven by the magnetic field gradient/curvature and pressure gradients. In the following analyses, we consider the Z pinch configuration as a prototype of magnetospheric plasmas where the magnetic field is approximated by a circle. The fast dynamics along the field line is ignored.

Figure 1.

A schematic view of the magnetospheric configuration governed by planetary dipolar magnetic fields. Low-frequency fluctuations in the equatorial plane are considered. The magnetic field gradient/curvature and pressure gradients drive drift waves. The electron (ion) drift is directed eastward (westward). The quasi-neutrality constraint prevents the entropy mode from travelling in a certain direction. However, the imbalance due to temperature difference allows the mode to travel and grow.

In the Z pinch, a plasma in a cylinder is confined by the azimuthal magnetic field |$\boldsymbol {B}=B_{\theta }(r)\nabla \theta$|⁠. We take a circular magnetic field of the Z pinch with the curvature radius given by |$R_{\mathrm{c}}$| and solve a local electrostatic gyrokinetic model for electrons and single-species ions. The background plasma is characterized by the Maxwellian distribution function,

$$\begin{eqnarray} f_{0s}=n_{0} \left(\frac{m_{s}}{2\pi T_{0s}}\right)^{3/2} \exp \left(-\frac{m_{s} v^2}{2T_{0s}}\right), \end{eqnarray}$$

(1)

where |$n_{0}$| is the density, |$m_{s}$| is the mass, and |$T_{0s}$| is the temperature. The subscript |$s={\mathrm{i}}, {\mathrm{e}}$| signifies the species. Since there is no variation along the field line, it is sufficient to consider a two-dimensional slab perpendicular to the magnetic field. Hereafter, we take the Cartesian coordinate where z is the direction along the magnetic field and x is the direction of the background gradient, |$\boldsymbol {B}=B_{0}(x) \nabla z$| and |$(\nabla B_{0})/B_{0} = - (\mathrm{d}\ln B_{0}/\mathrm{d}x)\nabla x$| at the local position |$x=R_{\mathrm{c}}$|⁠.

The gyrokinetic model (Frieman & Chen 1982; Howes et al. 2006) is a suitable model for describing low-frequency microturbulence under ambient magnetic fields. It governs the evolution of the non-Boltzmann part of the perturbed distribution function |$h_s$| where the total perturbation |$\delta f_s = - (q_s \phi /T_{0s}) f_{0s} + h_{s}$| as

$$\begin{eqnarray} &&\frac{\partial h_{s}}{\partial t} + (\boldsymbol {v}_{\mathrm{D},s}+\boldsymbol {v}_{\boldsymbol {E},s}) \cdot \nabla h_{s} \\ && =\left(\frac{\partial \langle \phi \rangle _{\boldsymbol {R}_{s}}}{\partial t} + \boldsymbol {v}_{\ast ,s}^{\mathrm{T}} \cdot \nabla \left\langle \phi \right\rangle _{\boldsymbol {R}_s} \right) \frac{q_{s} f_{0s}}{T_{0s}} + C(h_{s}), \end{eqnarray}$$

(2)

complemented by the quasi-neutrality condition determining the electrostatic potential |$\phi$|⁠,

$$\begin{eqnarray} \sum _{s} q_{s} \delta n_{s} = 0, \end{eqnarray}$$

(3)

where |$\delta n_{s}=\int \delta f_{s} \mathrm{d}\boldsymbol {v}$| is the density perturbation, |$q_{s}$| is the charge, |$\boldsymbol {v}_{\mathrm{D},s}$|⁠, |$\boldsymbol {v}_{\ast ,s}^{\mathrm{T}}$| and |$\boldsymbol {v}_{\boldsymbol {E},s}=-\nabla \langle \phi \rangle _{\boldsymbol {R}_{s}}/B_{0}\times \nabla z$| are, respectively, the magnetic, diamagnetic, and |$\boldsymbol {E}\times \boldsymbol {B}$| drift velocities. The angle bracket denotes the gyro-averaging operation. We assume the collision term C, including an artificial hyperdiffusion (Ricci, Rogers & Dorland 2006a), is small. It is ignored in most of the following analyses but employed only in non-linear simulations just to achieve steady states. We decompose variables in Fourier series, |$h_{s}(\boldsymbol {R}_{s},v_{\parallel },v_{\perp },t) = \sum _{\boldsymbol {k}} h_{\boldsymbol {k},s}(v_{\parallel },v_{\perp }) e^{-\mathrm{i}(\omega (\boldsymbol {k}) t-\boldsymbol {k}\cdot \boldsymbol {R}_{s})}$|⁠, |$\phi (\boldsymbol {r},t)= \sum _{\boldsymbol {k}} \phi _{\boldsymbol {k}} e^{-\mathrm{i}(\omega (\boldsymbol {k}) t-\boldsymbol {k}\cdot \boldsymbol {r})}$|⁠, where |$\boldsymbol {r}=\boldsymbol {R}_{s}-\boldsymbol {v}\times \nabla z/\Omega _{\mathrm{c},s}$| with |$\boldsymbol {r}$|⁠, |$\boldsymbol {R}_{s}$| being the particle and gyrocentre coordinates. The collisionless linear solution is given by

$$\begin{eqnarray} h_{\boldsymbol {k},s} = \frac{\omega -\omega _{\ast ,s}^{\mathrm{T}}}{\omega -\omega _{\mathrm{D},s}} J_{0}(\alpha _{s}) \frac{q_{s} \phi _{\boldsymbol {k}}}{T_{0s}} f_{0s}. \end{eqnarray}$$

(4)

The drift frequencies are defined by

$$\begin{eqnarray} \omega _{\ast ,s}^{\mathrm{T}} = \boldsymbol {k}\cdot \boldsymbol {v}_{\ast ,s}^{\mathrm{T}} = \omega _{\ast ,n,s} (1 + \eta _{s} (v^2/v_{\mathrm{th},s}^2-3/2)), \end{eqnarray}$$

(5)

$$\begin{eqnarray} \omega _{\mathrm{D},s} = \boldsymbol {k} \cdot \boldsymbol {v}_{\mathrm{D},s} = \omega _{\kappa ,s} v_{\parallel }^{2}/v_{\mathrm{th},s}^{2} + \omega _{\nabla B,s} v_{\perp }^{2}/(2 v_{\mathrm{th},s}^{2}), \end{eqnarray}$$

(6)

where the velocity-independent drift frequencies are, respectively, given by

$$\begin{eqnarray} \omega _{\ast ,n,s} = -\frac{1}{2} k_{y} \rho _{s} v_{\mathrm{th},s}/L_{n}, \end{eqnarray}$$

(7)

$$\begin{eqnarray} \omega _{\kappa ,s} = -k_{y} \rho _{s} v_{\mathrm{th},s}/R_{\mathrm{c}}, \end{eqnarray}$$

(8)

$$\begin{eqnarray} \omega _{\nabla B,s} = -k_{y} \rho _{s} v_{\mathrm{th},s}/L_{B}. \end{eqnarray}$$

(9)

The Larmor radius |$\rho _{s}=v_{\mathrm{th},s}/\Omega _{\mathrm{c},s}$|⁠, the thermal velocity |$v_{\mathrm{th},s}=\sqrt{2T_{0s}/m_{s}}$|⁠, and the cyclotron frequency |$\Omega _{\mathrm{c},s}=q_{s} B_{0}/m_{s}$|⁠. Note that the cyclotron frequency is defined including the sign of charge. The electron drift frequencies are positive, indicating electrons drift in the positive y direction. The background density and magnetic field gradient scales are |$L_{n}\equiv -\mathrm{d}\ln n_{0}/\mathrm{d}x$| and |$L_{B}\equiv -\mathrm{d}\ln B_{0}/\mathrm{d}x$| and the magnetic field curvature radius is |$R_{\mathrm{c}}$|⁠. For the electrostatic case, the MHD equilibrium condition demands |$L_{B}=R_{\mathrm{c}}$|⁠. In this case, we write |$\omega _{\kappa ,s}=\omega _{\nabla B,s}\equiv \omega _{B,s}$|⁠. The ratio of temperature and density gradients is |$\eta _{s}=\mathrm{d}\ln T_{0s}/\mathrm{d}\ln n_{0}$|⁠. In the following analyses, we ignore |$\eta _{s}$| for simplicity. The argument of Bessel function |$J_{0}$| is |$\alpha _{s}=k v_{\perp }/\Omega _{\mathrm{c},s}$| with |$k=|\boldsymbol {k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$|⁠.

We introduce the standard notation of electrostatics in a dielectric medium. The polarization is defined by |$q_{s} \delta n_{s} = -\nabla \cdot \boldsymbol {P}_{s}$| where |$\boldsymbol {P}_{s}=\epsilon _{0} \chi _{s} \boldsymbol {E}$| using the electric susceptibility |$\chi _{s}$| (⁠|$\epsilon _{0}$| is the vacuum permittivity) with |$\boldsymbol {E}=-\nabla \phi$| being the electrostatic field. Since the polarization is time-dependent, it carries the polarization current |$\boldsymbol {j}_{\mathrm{pol},s}=\partial \boldsymbol {P}_{s}/\partial t = \sigma _{s} \boldsymbol {E}$|⁠. The associated conductivity is given by |$\sigma _{s} = - \mathrm{i}\omega \epsilon _{0} \chi _{s}$|⁠. Using the solution of the gyrokinetic equation (4) for the density perturbation, we obtain the susceptibility as

$$\begin{eqnarray} \chi _{s} = \frac{1}{k^{2} \lambda _{\mathrm{D},s}^{2}} \left(1-\frac{\omega _{\ast ,n,s}-\omega }{\omega _{B,s}} Q(\omega /\omega _{B,s},k \rho _{s})\right), \end{eqnarray}$$

(10)

where |$\lambda _{\mathrm{D},s}\equiv \sqrt{\epsilon _{0} T_{0s}/(q_{s}^{2}n_{0})}$| is the Debye length and the integral Q is given by

$$\begin{eqnarray} && Q(\omega /\omega _{B,s}, k \rho _s) \\ && = \frac{2}{\sqrt{\pi }} \int \frac{e^{-(\hat{v}_{\parallel }^2 + \hat{v}_{\perp }^{2})}}{\hat{v}_{\parallel }^2 + \hat{v}_{\perp }^2/2- \omega /\omega _{B,s}} J_{0}^{2} (k \rho _{s} \hat{v}_{\perp }) \hat{v}_{\perp } \mathrm{d}\hat{v}_{\perp } \mathrm{d}\hat{v}_{\parallel }. \end{eqnarray}$$

(11)

From the quasi-neutrality, which is the small |$k\lambda _{\mathrm{D},s}$| limit of Gauss’s law, the dispersion relation is given by |$\chi \equiv \chi _{\mathrm{e}} + \chi _{\mathrm{i}} = 0$| or |$\sigma \equiv \sigma _{\mathrm{e}} + \sigma _{\mathrm{i}} = 0$|⁠.

The integral Q can be analytically evaluated if |$k\rho _{s}\ll 1$| (Ricci et al. 2006b). If we ignore the FLR correction, Q is simply equal to |$-\Xi (\mathrm{i}\sqrt{-\omega /\omega _{B,s}})^2$| where |$\Xi (z)=\mathrm{i}\sqrt{\pi } e^{-z^2} \mathrm{erfc}(-\mathrm{i}z)$| is the plasma dispersion function, |$\mathrm{erfc}(z)$| is the complementary error function.

The total energy of the system per unit volume is given by the sum of the energy of each species, |$W=\sum _{s}W_{s}$|⁠, where

$$\begin{eqnarray} W_{s} = \frac{1}{V} \int \frac{T_{0s} \delta f_{s}^{2}}{2f_{0s}} \mathrm{d}\boldsymbol {v} \mathrm{d}\boldsymbol {r} = \sum _{\boldsymbol {k}} \int \frac{T_{0s}}{2f_{0s}} \left| \delta f_{\boldsymbol {k},s} \right|^2 \mathrm{d}\boldsymbol {v}, \end{eqnarray}$$

(12)

with |$V=L_{x}L_{y}$| and |$L_{x}$|⁠, |$L_{y}$| being the box sizes in the x, y directions. The time evolution of the particle energy is given by

$$\begin{eqnarray} \frac{\mathrm{d} W_{s}}{\mathrm{d} t} = H_{s} + I_{s}-D_{s}, \end{eqnarray}$$

(13)

$$\begin{eqnarray} H_{s} = \frac{1}{V}\int q_{s} h_{s} \frac{\partial \langle \phi \rangle _{\boldsymbol {R}_{s}}}{\partial t} \mathrm{d}\boldsymbol {v} \mathrm{d}\boldsymbol {r}, \end{eqnarray}$$

(14)

$$\begin{eqnarray} I_{s} = \frac{T_{0s}}{L_{n}} \varGamma _{s} = \frac{T_{0s}}{L_{n}} \frac{1}{V} \int h_{s} \boldsymbol {v}_{\boldsymbol {E},s}\cdot \nabla x \mathrm{d}\boldsymbol {v} \mathrm{d}\boldsymbol {r}, \end{eqnarray}$$

(15)

$$\begin{eqnarray} D_{s} = - \frac{1}{V} \int \frac{T_{0s} h_{s}}{f_{0s}} C(h_{s}) \mathrm{d}\boldsymbol {v} \mathrm{d}\boldsymbol {r}, \end{eqnarray}$$

(16)

where |$H_{s}$| is the energy exchange between species, |$I_{s}$| is the energy injection, |$\varGamma _{s}$| is the particle flux and |$D_{s}$| is the collisional dissipation (entropy production). We can show that |$\sum _{s} H_{s} = 0$|⁠. Therefore, in the non-linearly saturated state, |$\sum _{s} (I_{s}-D_{s}) = 0$|⁠. The energy exchange term is known as the turbulent heating as it contributes to the increase of background temperature evolution in a longer time-scale as (Sugama et al. 1996; Howes et al. 2006; Candy 2013; Barnes et al. 2018),

$$\begin{eqnarray} \frac{3}{2} n_{0} \frac{\mathrm{d}T_{0s}}{\mathrm{d}t^{\prime }} + \frac{\partial \overline{\mathcal {Q}}_{s}}{\partial x^{\prime }} = \overline{H_{s}}, \end{eqnarray}$$

(17)

where |$t^{\prime }$|⁠, |$x^{\prime }$| denote the macroscale coordinates compared with the microscopic gyrokinetic t, x, the overline denotes a long-time average in t and |${\mathcal {Q}}_{s}$| is the turbulent heat flux.

If we substitute the collisionless linear solution (4) to (12), we find

$$\begin{eqnarray} 2 \Im (\omega) W_{\boldsymbol {k},s} = \epsilon _{0}|\boldsymbol {E}_{\boldsymbol {k}}|^{2} \left(\Re \left(\frac{\sigma _{s}}{\epsilon _{0}} \right)-\omega _{\ast ,n,s} \Im (\chi _{s}) \right) \end{eqnarray}$$

(18)

satisfied by each |$\boldsymbol {k}$| where the first and second terms on the right-hand side, respectively, correspond to the energy exchange and injection. (In fact, the energy exchange is equivalent to the work done on the species s, |$\propto \sigma _{s}|\boldsymbol {E}|^{2}$|⁠). It is evident that the species having negative conductivity must extract the energy from the background to grow its energy (i.e. |$\Im (\omega)\gt 0$|⁠).

3 STABILITY ANALYSES IN TERMS OF ENERGETICS

We revisit the linear stability analyses of the entropy mode by Ricci et al. (2006a and b) and identify the free energy source for the instability according to the prescription given in Hasegawa (1975). By looking at the direction of decreasing the free energy, we predict the state to which the system tries to relax.

As in Ricci et al. (2006b), we neglect the electron FLR effect, |$k \rho _{\mathrm{e}} \rightarrow 0$|⁠, but keep the small but finite ion FLR effect, |$k \rho _{\mathrm{i}}\ll 1$|⁠. We also assume |$\omega /\omega _{B,s}\ll 1$| to reduce the plasma dispersion function, |$\Xi (z)\approx \mathrm{i}\sqrt{\pi }-2 z$|⁠. To consider the marginal state, we write |$\omega =\omega _{0} + \mathrm{i}\Im (\omega)$| with |$\Im (\omega)\ll \omega _{0}$|⁠. By Taylor expanding the dispersion relation |$\sigma = 0$| around |$\omega =\omega _{0}$|⁠, the solution must satisfy

$$\begin{eqnarray} \Im (\sigma (\omega =\omega _{0})) = - \omega _{0} \epsilon _{0} \Re (\chi (\omega =\omega _{0}))=0, \end{eqnarray}$$

(19)

$$\begin{eqnarray} \Im (\omega) = -\frac{\Re (\sigma (\omega =\omega _{0}))}{\left(-\partial \Im (\sigma)/\partial \omega |_{\omega =\omega _{0}}\right)}. \end{eqnarray}$$

(20)

The numerator and denominator in (20) indicate, respectively, the dissipation due to wave–particle interactions and the wave energy equals to |$-\left.\partial \Im (\sigma)/\partial \omega \right|_{\omega =\omega _{0}} \approx \epsilon _{0} \left.\partial \Re (\omega \chi)/\partial \omega \right|_{\omega =\omega _{0}}$| (Landau, Lifshitz & Pitaevskii 1984). Therefore, (20) states that the coupling of the negative (positive) dissipation with the positive (negative) energy wave leads to instability. This generalizes the statement that the negative energy is the necessary condition for instability in the systems without wave–particle interactions or other loss mechanisms.

Care must be taken when approximating the square roots to obtain the solution. According to the solution given in Ricci et al. (2006b), |$\omega _{0} \gt 0$| for |$\tau /Z\lt 1$| where |$\tau \equiv T_{0\mathrm{i}}/T_{0\mathrm{e}}$| is the temperature ratio and |$Z\equiv -q_{\mathrm{i}}/q_{\mathrm{e}}$| is the ion charge number (we always set |$Z=1$|⁠). For positive |$\Im (\omega)$|⁠, |$\sqrt{-\omega /\omega _{B,\mathrm{i}}}\approx \sqrt{(Z/\tau)(\omega _{0}/\omega _{B,\mathrm{e}})}$|⁠, and |$\sqrt{-\omega /\omega _{B,\mathrm{e}}}\approx -\mathrm{i}\sqrt{\omega _{0}/\omega _{B,\mathrm{e}}}$|⁠. Then, we obtain

$$\begin{eqnarray} \frac{\omega _{0}}{\omega _{\ast ,n,\mathrm{e}}} & =& \frac{\tau }{Z} \frac{1}{2\pi (R_{\mathrm{{c}}}/L_{n})^{3}} \\ && \times \left[\left(1+\frac{\tau }{Z}\right)\left(\pi \frac{R_{\mathrm{c}}}{2L_{n}}-1\right)- \frac{k^{2}\rho _{\mathrm{i}}^{2}}{2} \frac{R_{\mathrm{c}}}{2L_{n}}(\pi -2) \right]^{2}, \end{eqnarray}$$

(21)

$$\begin{eqnarray} \frac{\Im (\omega)}{\omega _{\ast ,n,\mathrm{e}}} = 2 \left(\frac{\tau }{Z}\right)^{3/2} \frac{\omega _{0}}{\omega _{\ast ,n,\mathrm{e}}}, \end{eqnarray}$$

(22)

which is the solution for |$\tau /Z\ll 1$| of Ricci et al. (2006b). A positive imaginary part results from (20) because

$$\begin{eqnarray} \Re (\sigma _{\mathrm{e}}(\omega =\omega _{0})) = - \omega _{0} \frac{\epsilon _{0}}{k^{2}\lambda _{\mathrm{D},\mathrm{e}}^{2}} 4 \sqrt{\pi } \frac{\omega _{\ast ,n,\mathrm{e}}}{\omega _{B,\mathrm{e}}} \sqrt{\frac{\omega _{0}}{\omega _{B,\mathrm{e}}}} \lt 0, \end{eqnarray}$$

(23)

$$\begin{eqnarray} -\left.\frac{\partial \Im (\sigma _{\mathrm{i}})}{\partial \omega }\right|_{\omega =\omega _{0}} = \frac{1}{2}\sqrt{\frac{Z}{\tau }} \frac{\epsilon _{0}}{k^{2}\lambda _{\mathrm{D},\mathrm{i}}^{2}} 4 \sqrt{\pi } \frac{\omega _{\ast ,n,\mathrm{i}}}{\omega _{B,\mathrm{i}}} \sqrt{\frac{\omega _{0}}{\omega _{B,\mathrm{e}}}} \gt 0. \end{eqnarray}$$

(24)

[Note that the electron negative dissipation, (23), also indicates that the electron energy is negative. Under the present assumption, |$\omega _{0}\gg \Im (\omega)\gt 0$|⁠, the negative dissipation and negative energy are equivalent.] Therefore, we find that, for |$\tau /Z\lt 1$|⁠, the wave propagates in the electron direction and the system becomes unstable due to the coupling of positive energy ions and negative dissipation electrons. By performing the same analysis for |$\tau /Z\gt 1$|⁠, we find that the positive energy electrons couple with the negative dissipation ions to lead to instability. According to the direction of wave propagation, the roles played by ions and electrons flip. Fig. 2 shows the conductivity |$\sigma _{s}$| for |$\tau =0.5, 1, 2$| and |$L_{n}/R_{\mathrm{c}}=0.5$| obtained from the linear gyrokinetic solutions including full FLR of both ions and electrons. We show |$\Re (\sigma _{s})/\epsilon _{0} (L_{n}/v_{\mathrm{th,i}}) (\lambda _{\mathrm{D},s}/\rho _{\mathrm{Se}})^{2}$| against the wavenumber |$k\rho _{\mathrm{Se}}$|⁠. For |$\tau \lt 1$|⁠, we confirm that the electron conductivity is negative as calculated above. For |$\tau \gt 1$|⁠, the ion conductivity is negative only for small |$k\rho _{\mathrm{Se}}$|⁠, but it becomes positive as |$k\rho _{\mathrm{Se}}$| gets large, maybe because of the stabilizing influence of the FLR effect. The range of negative ion conductivity broadens as |$L_{n}/R_{\mathrm{c}}$| becomes large. For |$\tau =1$|⁠, the eigenvalue |$\omega$| is pure imaginary, and the assumption that |$\omega _{0} \gg \Im (\omega)$| is violated. In this case, it is the electron conductivity that is negative as the numerical solution shows. We expect that the sign of conductivity changes slightly above |$\tau =1$|⁠.

Figure 2.

The conductivity calculated from linear gyrokinetic solutions. Depending on the temperature ratio, the conductivity of each species changes its sign. The negative conductivity species plays a destabilizing role. For |$\tau \ne 1$|⁠, the conductivity is large at small k.

4 NON-LINEAR SIMULATION OF TURBULENT HEATING

What are the non-linear consequences of these instability mechanisms? One may speculate that the sign of the conductivity corresponds to the heating direction, namely, for |$\tau /Z\lt 1$|⁠, the negative electron conductivity indicates the positive ion heating, and vice versa. If that is true, the plasma tends to thermally equilibrate between species to achieve an equal temperature state, i.e. |$\tau$| approaches to the unity. Note that this prediction is not trivial because the turbulent heating is non-linear and depends on the turbulent cascades in phase space. For example, if we consider a uniform background and ignore the effect of resonance along |$\omega =\omega _{\mathrm{D},s}$| (see 4), electrons and ions are energetically decoupled in the electrostatic plasmas (Schekochihin, Kawazura & Barnes 2019), thus non-linear energy exchange does not occur. The linear term opens a channel of energy exchange between species, pre-dominantly in the energy injection scale. Whether the energy exchange is significant or not depends on how broad the injection scale in wavenumber space is. To confirm the speculation that the entropy-mode turbulence drives the system toward the thermal equilibrium, we perform non-linear gyrokinetic simulations.

Numerical simulations are performed using the gs2 code (Barnes et al. 2024). We fix the ion charge |$Z=1$|⁠, mass ratio |$m_{\mathrm{i}}/m_{\mathrm{e}}=1836$|⁠, density gradient |$L_{n}/R_{\mathrm{c}}=0.5, 2/3$|⁠, and change the temperature ratio |$\tau$|⁠. We take those density gradient values as they are just above the MHD interchange instability threshold |$\approx 0.3$| (Ricci et al. 2006a and b). [From the observations in the Earth’s magnetosphere (e.g. Denton et al. 2004), |$L_{n}/R_{\mathrm{c}}\approx 1$| is a reasonable choice.] The collisionality for ions and electrons are both small, |$\nu L_{n}/v_{\mathrm{th,i}}=10^{-2}$|⁠. We include the hypervisicosity (Ricci et al. 2006a) to achieve a saturated state. The amplitude is |$D_{\mathrm{hyp}}L_{n}/v_{\mathrm{th,i}}=1$| and the cut-off wavenumber is |$\approx 0.8 k_{\mathrm{max}}$|⁠. (The maximum wavenumber in the simulation is |$k_{\mathrm{max}}\rho _{\mathrm{Se}} \approx 2\pi (\rho _{\mathrm{Se}}/L) (N/3)$| where L and N are the box size and grid number.) The simulation domain is |$L_{x}/\rho _{\mathrm{Se}}=L_{y}/\rho _{\mathrm{Se}}=40\sqrt{2}\pi$| and the numerical resolution is |$(N_{x},N_{y},N_{\lambda },N_{E})=(128,128,8,32)$| where |$\lambda \equiv v_{\perp }^{2}/v^{2}$| and |$E\equiv v^{2}/v_{\mathrm{th},s}^2$| are the velocity space coordinates. To measure the heating in simulations, we use the symmetric form proposed by Candy (2013),

$$\begin{eqnarray} H_{s} = \frac{1}{2V} \int q_{s}\left(h_{s} \frac{\partial \langle \phi \rangle _{\boldsymbol {R}_{s}}}{\partial t}-\frac{\partial h_{s}}{\partial t} \langle \phi \rangle _{\boldsymbol {R}_{s}}\right) \mathrm{d}\boldsymbol {v} \mathrm{d}\boldsymbol {r}, \end{eqnarray}$$

(25)

which is equivalent to (14) if the long-time average is taken. Fig. 3 shows the time-averaged ion heating |$\overline{H_{\mathrm{i}}}$| in the saturated states, normalized by |$(\rho _{\mathrm{Se}}/L_{n})^{2} n_{0} T_{0\mathrm{e}} v_{\mathrm{th,i}}/L_{n}$|⁠. We first discuss the strongly driven case by the steep density gradient, |$L_{n}/R_{\mathrm{c}}=0.5$|⁠. As expected, the ion heating is positive when |$\tau$| is small and decreases as |$\tau$| increases. The sign changes around |$\tau \approx 1$|⁠. The turbulent heating is remarkable only when turbulence is strongly driven. When the density gradient becomes shallower (⁠|$L_{n}/R_{\mathrm{c}}=2/3$|⁠), the zonal flows self-generated from turbulence weaken turbulent fluctuations and resultant heating. [We still see the positive (negative) ion heating for |$\tau \lesssim 1$| (⁠|$\tau \gtrsim 1$|⁠).] We have confirmed the amplitude of the heating rate is generally in the same order as the energy injection.

Figure 3.

The time-averaged ion heating |$\overline{H_{\mathrm{i}}}$| in the saturated states obtained from non-linear gyrokinetic simulations. For the strongly driven case, |$L_{n}/R_{\mathrm{c}}=0.5$|⁠, the ion heating is negative for |$\tau \gt 1$| and positive for |$\tau \lt 1$|⁠. Therefore, ions will be cooled (heated) for |$\tau \gt 1$| (⁠|$\tau \lt 1$|⁠) to approach the equal temperature state. For the weakly driven case, |$L_{n}/R_{\mathrm{c}}=2/3$|⁠, turbulence is suppressed by self-generated zonal flows and so as the turbulent heating.

Fig. 4 shows perturbed distribution functions of ions and electrons for |$L_{n}/R_{\mathrm{c}}=0.5$| and |$\tau =0.5, 2$|⁠, taken at randomly sampled spatial points in the non-linearly saturated states. The black solid line shows the line of |$\omega _{\mathrm{D},s}=\omega _{B,s}$| for reference, while the actual resonance occurs along |$\omega = \omega _{\mathrm{D},s}$|⁠. We observe the signature of phase mixing in both ions and electrons, namely, the oscillatory structures in both |$v_{\parallel }$| and |$v_{\perp }$| directions. The structures are fluctuating in space and time, therefore must be subject to more sophisticated analyses, such as the spectral decomposition (Tatsuno et al. 2009; Kawazura et al. 2019; Meyrand et al. 2019; Nastac et al. 2024) or the field–particle correlation analyses (Klein & Howes 2016; Howes, Klein & Li 2017) to determine their secular effects on the energy exchange and dissipation.

Figure 4.

Sample perturbed distribution functions of ions (left) and electrons (right) for |$L_{n}/R_{\mathrm{c}}=0.5$|⁠, |$\tau =0.5$| (top) and |$\tau =2$| (bottom) cases. The distribution functions are taken at randomly chosen spatial points in the non-linearly saturated states. At some point in space and time, the distribution functions show unique structures typical of phase mixing. In the bottom right panel, the distribution function is close to Maxwellian, while in the other three panels, the distribution functions are oscillating in both |$v_{\parallel }$| and |$v_{\perp }$| directions indicating phase mixing.

In the present system, the phase mixing due to the resonance at |$\omega =\omega _{\mathrm{D},s}$| causing the continuum damping (Kim et al. 1994; Kuroda et al. 1998) as well as the non-linear phase mixing due to the FLR effect (Dorland & Hammett 1993; Tatsuno et al. 2009) can occur. Those phase space cascades transfer the injected energy to small-scale fluctuations, which are eventually dissipated by collisions. The amount of dissipation crucially depends on the dominant cascade mechanisms, which is also related to the energy exchange as |$H_s = D_s-I_s$| in the steady state. The cascade in velocity space is intrinsically two-dimensional (both |$v_{\perp }$| and |$v_{\parallel }$| are involved). Reliable numerical simulations having a sufficient resolution in the phase space and without using artificial hyperviscosity are computationally challenging, which we do not go into detail here.

5 CONCLUSION

Turbulence is a symptom of violent energy relaxation toward an equilibrium state. An excess free energy characterizing deviation from the equilibrium state is rapidly released via turbulence. In other words, finding the instability source directly indicates what is the preferable state. We have considered the entropy mode turbulence in the Z pinch configuration as an example. It has been shown that the drift wave propagation direction and responsible species for destabilization change according to the temperature ratio |$\tau$|⁠. If the ion temperature is larger than that of electrons, ions take the energy from the background to heat electrons, and vice versa. Therefore, the plasma tries to equalize the temperature even without collisions. Non-linear simulations confirm that the preferable temperature ratio is |$\tau \approx 1$|⁠. Although the instability and resultant turbulent anomalous transport were well known, it has been newly discovered by examining the instability source that the rapid energy exchange process between species is also operating in magnetospheric turbulence.

We have examined heating in a strongly turbulent state driven by the steep density gradient. However, if the turbulence is weak, the zonal flows regulate transport and also heating. The zonal flow generation and turbulence suppression in the Z pinch were intensively studied in Hallenbert & Plunk (2022) to predict the precise condition for the onset of turbulence. To comprehensively understand the global self-organized state in the magnetosphere, the transport and mesoscale structure formation, as well as the newly proposed energy exchange mechanism, must be taken into account simultaneously.

The entropy-mode turbulence and energy exchange between species seem rather universal in the magnetospheric environment. Yet, how dominant this mechanism is over other phenomena must be examined via space and experimental observations. Here, we roughly estimate the relevance to the Earth’s magnetosphere. In the radiation belt having its scale of about |$10^6-10^7$| m, roughly the MeV order of protons and electrons are confined (see, e.g. Wolf 1995). For |$k\rho _{\mathrm{Se}}=O(1)$| fluctuations, the drift frequency is about |$1-10$| Hz, which is about the upper limit of the ultra-low-frequency (ULF) fluctuations. The observed magnetospheric plasmas are typically in high-|$\beta$| (where |$\beta$| is the ratio of thermal and magnetic pressures) and are associated with electromagnetic fluctuations in the ULF range (Kivelson 1995). As discussed in Simakov et al. (2002), the entropy mode is purely electrostatic and does not involve magnetic fluctuations. Therefore, the entropy mode may co-exist with electromagnetic fluctuations. Our preliminary calculations show that the entropy mode still exists in high-|$\beta$| plasmas. We leave more detailed studies on the coupling between the entropy mode and the magnetic fluctuations (Alfvén waves along the field lines) for future consideration.

Although it is not a target of this study, it is worth mentioning the application of the present mechanism to a certain type of accretion flows. [A possible application may be magnetically arrested disks (Narayan, Igumenshchev & Abramowicz 2003) where a magnetically dominated region is formed to disrupt flows directly accreting to the central object.] In hot accretion flows, plasmas are usually assumed to be two-temperature, namely, ions are significantly hotter than electrons, because no known mechanism exists that transfers energy from ions to electrons. Only the possible mechanism proposed by Begelman & Chiueh (1988) due to turbulent heating may not be sufficient to eliminate the two-temperature nature (Narayan & Yi 1995). The present mechanism is similar to that of Begelman & Chiueh (1988) since the drift-type instability is the source of turbulence in both of mechanisms. For comparison, we derive the electron heating rate for the |$\tau /Z\gt 1$| case using the mixing length argument that estimates the saturated amplitude of turbulence. By equating the linear growth rate with the non-linear |$\boldsymbol {E}\times \boldsymbol {B}$| time-scale, we obtain the saturated amplitude of the electrostatic potential as |$\phi /B_{0}\sim \Im (\omega)/k^2$| with |$\Im (\omega) \sim \omega _{\ast ,n,\mathrm{i}}$| and |$k \sim \rho _{\mathrm{i}}^{-1}$|⁠. Then, we obtain the heating rate of electrons (equals to the cooling rate of ions),

$$\begin{eqnarray} H_{\mathrm{e}} = - H_{\mathrm{i}} \sim n_{0} T_{0\mathrm{i}}|\omega _{\ast ,n,\mathrm{i}}|^{3} \Omega _{\mathrm{c,i}}^{-2} \end{eqnarray}$$

(26)

which is similar to that of Begelman & Chiueh (1988), but the power indices are slightly modified. Note that the entropy mode is only active when the temperature difference is not so large and is stabilized for |$\tau /Z\gg 1$|⁠. There might be a critical value of |$\tau /Z$| below which the two-temperature assumption breaks. To determine whether the present mechanism works or not, specific conditions of plasmas must be carefully considered.

ACKNOWLEDGEMENTS

The author would like to thank W. Dorland for his encouraging comments. This work was supported by Japan Society for Promotion of Science (JSPS) KAKENHI Grant Number JP22K03568. Numerical simulations were performed on ‘Plasma Simulator’ (NEC SX-Aurora TSUBASA) of National Institute for Fusion Science (NIFS) with the support and under the auspices of the NIFS Collaboration Research programe (NIFS22KISS019).

DATA AVAILABILITY

The data that support the findings of this study are openly available in Zenodo at 10.5281/zenodo.12663811 (Numata 2024).