1. Introduction
A diversity of astrophysical phenomena involves hot, diluted, and turbulent plasmas where ambient conditions make Coulomb collisions between particle species scarce, and therefore are weakly collisional or collisionless. This absence of collisions makes these systems prone to departures from thermodynamic equilibrium, giving rise to a plethora of kinetic scale phenomena that mediate the momentum and energy transport between scales. In this sense, kinetic scale processes become key to understanding the global thermodynamic evolution of these astrophysical settings. Some examples of astrophysical scenarios where these plasmas are present are the solar wind, low luminosity accretion disks around supermassive black holes (SMBH, e.g. Sgr A
$^*$
at the centre of our Galaxy, and the SMBH at the centre of M87) and the intracluster medium (ICM) of galaxy clusters. In addition to these scenarios, there are systems where these collisionless plasmas become relativistically hot. Some examples of this type include pulsar wind nebulae, relativistic jets from active galactic nuclei (e.g. blazars, radio-loud quasars) and cosmic rays.
An important feature of magnetised, collisionless systems is the development of a pressure anisotropy
$\Delta P_{\alpha } \equiv P_{\perp ,\alpha }-P_{\parallel ,\alpha }$
, where
$P_{\perp ,\alpha }$
,
$P_{\parallel ,\alpha }$
are the components of the pressure tensor of species
$\alpha$
perpendicular and parallel to the ambient magnetic field. This anisotropy in the pressure arises from the conservation of the particle adiabatic invariants
$\mathcal{M} = p_{\perp ,\alpha }^2/2B$
, d the so-called magnetic moment, where
$B$
denotes the magnetic field strength, and
$\mathcal{J} = \oint p_{\parallel ,\alpha }{\rm d}\ell$
, the longitudinal action, where
$p_{\perp ,\alpha }$
,
$p_{\parallel ,\alpha }$
are the perpendicular and parallel components of the momentum of species
$\alpha$
,
$\mathrm{d}\ell$
is the path differential, and the integral of the parallel action is taken between points of concentration of magnetic field lines, or so called bounce points. Written in this form, these invariants are conserved in both non-relativistic and relativistic regimes. The evolution of
$P_{\perp ,\alpha }$
and
$P_{\parallel ,\alpha }$
then decouples, and each component evolves independently.
The first description for the adiabatic evolution of the perpendicular and parallel pressures in a collisionless plasma was given by the influential paper by Chew, Goldberger & Low (Reference Chew, Goldberger and Low1956), now known as the double-adiabatic or CGL equations. These relations have been proven useful in a variety of contexts, e.g. for constructing fluid models with anisotropic closures to magnetohydrodynamic (MHD) equations (e.g. Hunana et al. Reference Hunana, Tenerani, Zank, Khomenko, Goldstein, Webb, Cally, Collados, Velli and Adhikari2019; Majeski, Kunz & Squire Reference Majeski, Kunz and Squire2024), as well as for contributing to our understanding of the temperature evolution in the solar wind (Matteini et al. Reference Matteini, Landi, Hellinger, Pantellini, Maksimovic, Velli, Goldstein and Marsch2007, Reference Matteini, Hellinger, Goldstein, Landi, Velli and Neugebauer2013). On the other hand, the validity of the CGL model is delimited by the validity of the conservation of adiabatic invariants. There exists a variety of kinetic processes that can break the adiabatic invariance, e.g. Coulomb collisions, heat fluxes and pressure anisotropy-driven microinstabilities. This way, CGL models have been one important first step into the richness of kinetic plasma processes.
Despite their great success, the CGL equations are only valid in the non-relativistic regime. Relativistic extensions to the double-adiabatic equations were proposed in Newcomb (Reference Newcomb1982), Holm & Kupershmidt (Reference Holm and Kupershmidt1986), Gedalin (Reference Gedalin1991) and Gedalin & Oiberman (Reference Gedalin and Oiberman1995). In Gedalin (Reference Gedalin1991) and Gedalin & Oiberman (Reference Gedalin and Oiberman1995), the most general state equations for an anisotropic collisionless plasmas were derived. In these equations,
$P_{\perp }$
and
$P_{\parallel }$
are expressed as a function of the density, magnetic field strength and the internal energy of the plasma
$\varepsilon$
, the latter calculated as the appropriate moment of a general, gyrotropic distribution function
$f(p_{\perp }^2,p_{\parallel }^2)$
, solution of the Vlasov equation for a well magnetised plasma, although no specific functional form for
$f(p_{\perp }^2,p_{\parallel }^2)$
was provided.
In this work, we directly solve the drift kinetic equation analytically, and obtain a general solution for a time-dependent, gyrotropic distribution function
$f(t,p_{\perp },p_{\parallel })$
, for a given initial, well-behaved distribution
$f_0(p_{\perp },p_{\parallel })$
. We then consider three different cases as initial conditions: non-relativistic Maxwellian, and relativistic and ultrarelativistic Maxwell–Jüttner thermal distributions, and calculate moments to obtain evolution equations for
$P_{\perp }$
and
$P_{\parallel }$
in each of the three different regimes. We directly recover the CGL equations for the non-relativistic case, and obtain new double-adiabatic evolution for the relativistic and ultrarelativistic cases, dependent on the density
$n(t)$
and field strength
$B(t)$
only. We confirm that the obtained double-adiabatic equations satisfy the general state equations proposed in Gedalin (Reference Gedalin1991) and Gedalin & Oiberman (Reference Gedalin and Oiberman1995) in both relativistic and ultrarelativistic regimes. We perform particle-in-cell (PIC) simulations with shear and compressive drivings to confirm our analytical results, and obtain remarkable agreement in both relativistic and ultrarelativistic regimes. The time-dependent distribution function obtained in the relativistic regime also constitutes a natural extension to an anisotropic Maxwell–Jüttner distribution, derived directly as a solution of the Vlasov equation. We confirmed that this anisotropic Maxwell–Jüttner distribution agrees with PIC simulations notably well.
A similar method of solving the Vlasov equation was presented in Zhdankin, Kunz & Uzdensky (Reference Zhdankin, Kunz and Uzdensky2023) for the case of an ultrarelativistic plasma with synchrotron cooling electrons. Incidentally, Zhdankin et al. (Reference Zhdankin, Kunz and Uzdensky2023) also presented a type of anisotropic distribution for ions, and obtained approximate relations for
$P_{\perp }$
and
$P_{\parallel }$
for a weak magnetic field strength. In this work we consider any level of anisotropy and magnetic field strength, allowing the results to be completely general.
This paper is organised as follows. In § 2 we present the main results of this work: we analytically solve the drift kinetic equation, and obtain time-dependent solutions for non-relativistic, relativistic and ultrarelativistic initial distributions. We then calculate moments of these solutions and obtain double-adiabatic evolution equations for
$P_{\perp }$
and
$P_{\parallel }$
in each regime. In § 3 we describe our simulation methods and set-up for our shearing and compressing PIC simulations. In § 4 we compare the analytical results of § 2 with PIC simulations of both shear and compressive drivings, for mildly relativistic and ultrarelativistic initial temperatures. We discuss possible future directions that this work can open in § 5. In § 6 we summarise our results and present our conclusions.
2. Theoretical basis
Consider a periodic, uniform, magnetised, collisionless flux tube of plasma of radius
$r$
and length
$l$
with no scattering mechanism active and zero heat fluxes. The evolution of the particle distribution function
$f_{\alpha }(t,p_{\perp },p_{\parallel })$
, which we assume gyrotropic, is given by the drift kinetic equation (Kulsrud Reference Kulsrud1983, Reference Kulsrud2005)
\begin{align} \frac {\partial f_{\alpha }}{\partial t} + \frac {{\rm d}p_{\perp }}{{\rm d}t}\frac {\partial f_{\alpha }}{\partial p_{\perp }} + \frac {{\rm d}p_{\parallel }}{{\rm d}t}\frac {\partial f_{\alpha }}{\partial p_{\parallel }} = 0, \end{align}
where
$p_{\perp }$
,
$p_{\parallel }$
are the perpendicular and parallel momenta with respect to the direction of the mean magnetic field, and
$\alpha$
denotes the particle species. We will assume that the density
$n_{\alpha }$
and magnetic field strength
$B$
vary very slowly in both spatial and temporal scales:
$L_n s_n, L_B s_B \gg v_{th,\alpha }$
, where
$L_n \equiv n/\boldsymbol{\nabla }n$
,
$L_B \equiv B/\boldsymbol{\nabla }B$
,
$s_n \equiv d\ln n_{\alpha }/{\rm d}t$
,
$s_B \equiv d\ln B/{\rm d}t$
and
$v_{th,\alpha }$
is the thermal velocity of the species
$\alpha$
. This way, particle free streaming cannot transport away variations of
$n$
and
$B$
, and we can neglect the spatial derivatives of
$f$
. The long parallel transit times implied under this assumption allow us to also neglect heat fluxes. We note that a complete, self-consistent treatment will have to include the convective derivative of
$f_{\alpha }$
. Particle drifts that can lead to finite off-diagonal terms in the pressure tensor and heat fluxes can, in general, exist, but consideration of them is beyond the scope of the present work.
Assuming that the spatial scales are much larger than the particle’s gyroradius,
$n/\boldsymbol{\nabla }n$
,
$B/\boldsymbol{\nabla }B$
$\gg R_{L,\alpha }$
and the plasma is well magnetised,
$\varOmega _{\alpha }/s \gg 1$
, where
$\varOmega _{\alpha } = eB/m_{\alpha }c^2$
is the non-relativistic cyclotron frequency of species
$\alpha$
, the following quantities become adiabatic invariants of the system: the particle magnetic moment
$\mathcal{M} = p_{\perp }^2/2B$
and the longitudinal action
$\mathcal{J} = \oint p_{\parallel } {\rm d}\ell$
, where the integral is taken between periodic points within the flux tube, and it becomes an invariant provided that the variation of the orbit is not substantial during a particle transit through the flux tube.Footnote
1
We can then use the conservation of magnetic moment, longitudinal action, particle number and magnetic flux to replace the time derivatives of
$p_{\perp },p_{\parallel }$
above, and obtain (see e.g. Lichko et al. Reference Lichko, Egedal, Daughton and Kasper2017)
\begin{align} \frac {\partial f_{\alpha }}{\partial t} + \frac {\dot {B}}{2B}p_{\perp }\frac {\partial f_{\alpha }}{\partial p_{\perp }} + \left ( \frac {\dot {n}_\alpha }{n_{\alpha }} - \frac {\dot {B}}{B} \right ) p_{\parallel }\frac {\partial f_{\alpha }}{\partial p_{\parallel }} = 0. \end{align}
Equation (2.2) is valid in both non-relativistic and relativistic regimes. If we go to spherical coordinates in momentum,
$p = (p_{\perp }^2 + p_{\parallel }^2)^{1/2}$
and
$ \mu = p_{\parallel }/p$
, the drift kinetic equation reads
\begin{align} &\frac {\partial f_{\alpha }}{\partial t} + \left ( \frac {\dot {n}_{\alpha }}{3n_{\alpha }}p + \!\left [ \frac {2}{3} \frac {\dot {n}_{\alpha }}{n_{\alpha }} - \frac {\dot {B}}{B}\right ]\frac {1}{2}(3\mu ^2 - 1)p\right )\frac {\partial f_{\alpha }}{\partial p} + \frac {3}{2}\left [ \frac {2}{3} \frac {\dot {n}_{\alpha }}{n_{\alpha }} - \frac {\dot {B}}{B}\right ]\mu (1 - \mu ^2)\frac {\partial f_{\alpha }}{\partial \mu } = 0.\nonumber \\& \end{align}
We remind the reader that our choice of variable
$\mu$
corresponds to the cosine of the pitch angle of the particle, not to be confused with the magnetic moment
$\mathcal{M}$
. Note that the term in square brackets can be written as
$(2/3)\dot {n}_{\alpha }/n_{\alpha } - \dot {B}/{B} = -\hat {b}\hat {b}:\boldsymbol{\nabla }\boldsymbol{u}_{\alpha } + (1/3)\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\alpha }$
, and
$(1/3)\dot {n}_{\alpha }/n_{\alpha } = -(1/3)\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\alpha }$
, so this combination of
$n$
and
$B$
variations effectively represents any external motion the plasma is subjected to. For convenience, define
$R \equiv n_{\alpha }^{-2/3}B$
, then we can simplify (2.3) to get
\begin{align} \frac {\partial f_{\alpha }}{\partial t} + \left ( \frac {\dot {n}_{\alpha }}{3n_{\alpha }}p - \frac {\dot {R}}{R}p\frac {1}{2}(3\mu ^2 - 1) \right )\frac {\partial f_{\alpha }}{\partial p} - \frac {3}{2}\frac {\dot {R}}{R}\mu (1 - \mu ^2)\frac {\partial f_{\alpha }}{\partial \mu } = 0. \end{align}
Given an initial distribution
$f_{\alpha ,0}(p,\mu )$
at some initial time
$t=t_0$
, (2.4) can be promptly solved by the method of characteristics to obtain a solution for
$f_{\alpha }(t,p,\mu )$
. Recognising that
$f_{\alpha }(t,p,\mu ) = f_{\alpha ,0}(p_0,\mu _0)$
and
\begin{align} p_0 &= p\left (\frac {n_{\alpha }}{n_{\alpha ,0}}\right )^{-1/3}\left \{ \left (\frac {R}{R_0} \right )^{-1} + \mu ^2\left [ \left (\frac {R}{R_0}\right )^2 - \left (\frac {R}{R_0}\right )^{-1}\right ] \right \}^{1/2}\nonumber \\[4pt] \mu _0 &= \left [ 1 + \frac {1-\mu ^2}{\mu ^2}\left (\frac {R}{R_0}\right )^{-3} \right ]^{-1/2}, \end{align}
we obtain
\begin{align} & f_{\alpha }(t,p,\mu ) = \nonumber\\& f_{\alpha ,0}\!\left ( p\left (\frac {n_{\alpha }}{n_{\alpha ,0}}\right )^{-1/3}\!\left \{ \left (\frac {R}{R_0} \right )^{-1}\! + \mu ^2\left [ \left (\frac {R}{R_0}\right )^2\! - \left (\frac {R}{R_0}\right )^{-1}\right ] \right \}^{1/2},\left [ 1\! + \frac {1-\mu ^2}{\mu ^2}\left (\frac {R}{R_0}\right )^{-3} \right ]^{-1/2} \right )\!, \end{align}
where
$R_0 = n_{\alpha ,0}^{-2/3}B_0$
are the initial values of
$n_{\alpha }(t)$
,
$B(t)$
at
$t=t_0$
. Given an initial distribution function
$f_{\alpha ,0}(p,\mu )$
, (2.6) then constitutes a completely general, time-dependent, homogeneous distribution function evolving in response to any specific
$B(t)$
,
$n_{\alpha }(t)$
variation. In this work, we will consider three different cases for
$f_{\alpha ,0}(p,\mu )$
: non-relativistic Maxwell–Boltzmann, relativistic Maxwell–Jüttner and ultrarelativistic Maxwell–Jüttner distributions.
One important feature of solution (2.6) will be the development of a pressure anisotropy
$P_{\perp } \neq P_{\parallel }$
, where
$P_{\perp }$
,
$P_{\parallel }$
are the pressures perpendicular and parallel to the direction of the mean magnetic field
$\boldsymbol{B}$
. We will see that, depending on the choice of
$f_{\alpha ,0}(p,\mu )$
, we can obtain evolution equations for
$P_{\perp }(t)$
and
$P_{\parallel }(t)$
as a function of
$B(t)$
and
$n(t)$
. These evolution equations will recover the classic CGL double-adiabatic equations (Chew et al. Reference Chew, Goldberger and Low1956) only in the non-relativistic regime, i.e. when considering a Maxwellian distribution for
$f_{\alpha ,0}(p,\mu )$
, and will depart from CGL when going to the relativistic and ultrarelativistic regimes (e.g. when considering a Maxwell–Jüttner distribution). Relevant additional studies using the solution distribution function (2.6), such as a linear stability analysis and different initial
$f_{\alpha ,0}(p,\mu )$
, will be deferred to future studies.
2.1. Non-relativistic case: Maxwell–Boltzmann distribution
Let us benchmark our results by considering first the case when the initial distribution
$f_{0}$
is a Maxwellian distribution function
\begin{align} f_{0}(v) = \frac {n_0}{\pi ^{3/2}v_{th}^3}\exp \!\left (-\frac {v^2}{v_{th}^2}\right )\!, \end{align}
where
$v_{th} = \sqrt {2k_BT_0/m}$
is the thermal velocity and
$T_{0}$
is the temperature, and we passed from momentum to velocity coordinates owing to the non-relativistic regime. For simplicity, let us drop the index
$\alpha$
, but bear in mind that this analysis is valid for any species. Then, (2.6) becomes
\begin{align} & f(t,v,\mu ) = \frac {n_0}{\pi ^{3/2}v_{th}^3}\exp \!\left (-\frac {v^2}{v_{th}^2} \left (\frac {n}{n_0}\right )^{-2/3}\left \{ \left (\frac {R}{R_0}\right )^{-1} + \mu ^2\left [ \left (\frac {R}{R_0}\right )^2 - \left (\frac {R}{R_0}\right )^{-1} \right ] \right \} \right )\!,\nonumber \\& \end{align}
where
$R/R_0 = (n/n_0)^{-2/3}(B/B_0)$
. Notably, we can do a coordinate transformation from
$(v,\mu )$
to
$(v_{\perp },v_{\parallel })$
and recover the well-known bi-Maxwellian distribution
\begin{align} f(t,v,\mu ) &= \frac {n_0}{\pi ^{3/2}v_{th}^3}\exp \!\left ( -\frac {v_{\perp }^2}{v_{th}^2(B/B_0)} - \frac {v_{\parallel }^2}{v_{th}^2(n/n_0)^{2}(B/B_0)^{-2}} \right )\!,\\[-10pt]\nonumber\end{align}
\begin{align} f(t,v,\mu ) &= \frac {n_0}{\pi ^{3/2}v_{th}^3}\exp \!\left ( -\frac {mv_{\perp }^2}{k_BT_{\perp }(t)} - \frac {mv_{\parallel }^2}{k_BT_{\parallel }(t)} \right )\!, \end{align}
where in the last step we replaced
$T_{\perp } = T_0(B/B_0)$
and
$T_{\parallel } = T_0(n/n_0)^{2}(B/B_0)^{-2}$
, as we naturally recover the expected, non-relativistic CGL evolution for
$T_{\perp }(t)$
and
$T_{\parallel }(t)$
(Chew et al. Reference Chew, Goldberger and Low1956). Indeed, if we take the
$P_{\perp }$
and
$P_{\parallel }$
non-relativistic moments of the distribution function (2.8), we also obtain the same result
\begin{align} P_{\perp } &= m\int \frac {v^2(1 - \mu ^2)}{2} f(t,v,\mu ) d^3 v = P_{0} \left (\frac {n}{n_0}\right )\left (\frac {B}{B_0}\right )\!,\\[-10pt]\nonumber\end{align}
\begin{align} P_{\parallel } &= m\int v^2\mu ^2 f(t,v,\mu ) d^3 v = P_{0} \left (\frac {n}{n_0}\right )^{3}\left (\frac {B}{B_0}\right )^{-2}, \end{align}
where
$P_0 = n_0k_BT_0$
, and the same can be obtained for the density
$n(t)$
. Given that our distribution functions are gyrotropic and neglecting finite Larmor radius effects, we can describe the pressure tensor by perpendicular and parallel scalars,
$P_{\perp }$
and
$P_{\parallel }$
. Then, we have recovered all the well-known results in the non-relativistic regime.Footnote
2
We will now see that when considering a Maxwell–Jüttner distribution, we obtain departures from the classic, non-relativistic CGL, and these new relations constitute the appropriate extension of the double-adiabatic equations to the relativistic and ultrarelativistic regimes for the initial distribution considered.
2.2. Ultrarelativistic case: Maxwell–Jüttner distribution
Let us start our relativistic extension of the double-adiabatic equations with the ultrarelativistic case
$\gamma = \sqrt {1 + (p/mc)^2} \gg 1$
, as it is analytically more tractable.
Consider now the case when the initial distribution function
$f_{0}$
is an ultrarelativistic Maxwell–Jüttner distribution
\begin{align} f_{0}(p) = \frac {n_0}{8\pi p_T^3}\exp \!\left (-\frac {p}{p_T}\right )\!, \end{align}
where
$p_T \equiv k_BT_0/c$
is the thermal momentum. Then, the distribution function (2.6) becomes
\begin{align} &f(t,p,\mu ) = \frac {n_0}{8\pi p_T^3}\exp \!\left ( -\frac {p}{p_T}\left (\frac {n}{n_0}\right )^{-1/3}\left (\frac {R}{R_0}\right )^{-1/2}\left \{ 1 + \mu ^2\left [ \left (\frac {R}{R_0}\right )^3 - 1 \right ] \right \}^{1/2} \right )\!,\nonumber \\&\end{align}
where
$R(t) \equiv n(t)^{-2/3}B(t)$
. The evolution of the perpendicular and parallel components of the pressure tensor are obtained by taking the corresponding moments of the distribution (2.15). Recalling that (e.g. Rezzolla & Zanotti Reference Rezzolla and Zanotti2013)
\begin{align} n &= \int f(t,p,\mu ) d^3p, \end{align}
and
\begin{align} T^{\mu \nu } &= c\int p^{\mu }p^{\nu } f \frac {d^3p}{p^0} \rightarrow c\int p^{\mu }p^{\nu } f \frac {d^3p}{p}, \end{align}
where
$p^{\mu }$
,
$p^{\nu }$
are the
$\mu$
and
$\nu$
components of the four-momentum, and the last expression in (2.16) is obtained by taking the ultrarelativistic limit
$p^0 = \gamma mc \rightarrow p$
. We then obtain, for the perpendicular pressure
\begin{align} P_{\perp } &= c\int \frac {p^2(1 - \mu ^2)}{2} f(t,p,\mu ) \frac {d^3 p}{p},\qquad\qquad\qquad\qquad\qquad\qquad \end{align}
\begin{align} P_{\perp }\! &= \frac {3}{4}P_0\! \left (\frac {n}{n_0}\right )^{4/3}\!\!\left (\frac {R}{R_0}\right )^2\! \left \{ \frac {(R/R_0)^3-2}{[(R/R_0)^3-1]^{3/2}}\arctan\! \left (\sqrt {(R/R_0)^3-1}\right ) + \frac {1}{(R/R_0)^3-1} \right \}\! ,\end{align}
for
$R\gt 1$
, and
\begin{align} P_{\perp }\! &= \frac {3}{4}P_0\!\left (\frac {n}{n_0}\right )^{4/3}\!\!\left (\frac {R}{R_0}\right )^2\! \left \{ \frac {2 - (R/R_0)^3}{[1-(R/R_0)^3]^{3/2}}\text{arctanh}\!\left (\sqrt {1 - (R/R_0)^3}\right ) - \frac {1}{1 - (R/R_0)^3} \!\right \}\!, \end{align}
for
$R\lt 1$
. For the parallel pressure
\begin{align} P_{\parallel } &= c\int p^2\mu ^2 f \frac {d^3 p}{p} ,\end{align}
\begin{align} P_{\parallel } &= \frac {3}{2}P_0 \left (\frac {n}{n_0}\right )^{4/3}\left (\frac {R}{R_0}\right )^2\left \{ \frac {\arctan \left ( \sqrt {(R/R_0)^3 - 1} \right )}{[(R/R_0)^3-1]^{3/2}} - \frac {1}{(R/R_0)^3[(R/R_0)^3-1]}\right \}\!,\nonumber \\& \end{align}
for
$R\gt 1$
, and
\begin{align} P_{\parallel } &= \frac {3}{2}P_0\left (\frac {n}{n_0}\right )^{4/3}\left (\frac {R}{R_0}\right )^2\left \{\frac {1}{(R/R_0)^3[1-(R/R_0)^3]} - \frac {\text{arctanh}\left (\sqrt {1 - (R/R_0)^3}\right )}{[1 - (R/R_0)^3]^{3/2}} \right \}\!,\nonumber \\& \end{align}
for
$R\lt 1$
. It can be shown that when
$B\rightarrow B_0$
and
$n\rightarrow n_0$
,
$P_{\perp } \rightarrow P_0$
and
$P_{\parallel } \rightarrow P_0$
, effectively recovering the initial isotropic state. An analogous calculation can be done to obtain the expected density evolution.
We can also calculate the internal energy density
$U$
in order to compare with the general state equations from Gedalin (Reference Gedalin1991)
\begin{align} U &= c\int p^2 f(t,p,\mu ) \frac {d^3p}{p} ,\end{align}
\begin{align} U &= \frac {3}{2}P_0\left (\frac {n}{n_0}\right )^{4/3}\left (\frac {R}{R_0}\right )^2\left \{\frac {\arctan \left (\sqrt {(R/R_0)^3-1}\right )}{\sqrt {(R/R_0)^3-1}} + \frac {1}{(R/R_0)^3}\right \} ,\end{align}
for
$R\gt 1$
, and
\begin{align} U &= c\int p^2 f(t,p,\mu ) \frac {d^3p}{p} ,\end{align}
\begin{align} U &= \frac {3}{2}P_0\left (\frac {n}{n_0}\right )^{4/3}\left (\frac {R}{R_0}\right )^2\left \{ \frac {\text{arctanh}\left (\sqrt {1 - (R/R_0)^3}\right )}{\sqrt {1 - (R/R_0)^3}} + \frac {1}{(R/R_0)^3} \right \} ,\end{align}
for
$R\lt 1$
(the two limits smoothly connect at
$R=1$
giving the correct value in the isotropic case). After a bit of algebra, it can be shown that
$P_{\perp }$
and
$P_{\parallel }$
satisfy the general state equations from Gedalin (Reference Gedalin1991)
\begin{align} P_{\parallel } &= n^2\frac {\partial }{\partial n}\left ( \frac {U}{n} \right )\!,\end{align}
\begin{align} P_{\perp } &= P_{\parallel } + nB\frac {\partial }{\partial B}\left (\frac {U}{n}\right )\!. \end{align}
Equations (2.18)–(2.22) constitute the extension of the double-adiabatic equations to the ultrarelativistic regime. We will see in § 4 that these relations correctly describe the evolution of
$P_{\perp }$
and
$P_{\parallel }$
in PIC simulations of shearing and compressing boxes.
2.3. Relativistic case: Maxwell–Jüttner distribution
Let us finally consider the full Maxwell–Jüttner distribution, valid for fully non-relativistic, moderately relativistic and ultrarelativistic temperatures
\begin{align} f_0(p) = \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p^2}{m^2c^2}}\right )\!, \end{align}
where
$\theta \equiv k_BT_0/mc^2$
is the normalised temperature and
$K_2(x)$
is the modified Bessel function of the second kind. Equation (2.6) then becomes
\begin{align} f(t,p,\mu ) &= \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\nonumber \\&\quad \exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p^2}{m^2c^2}\left (\frac {n}{n_0}\right )^{-2/3}\left (\frac {R}{R_0}\right )^{-1}\left \{ 1 + \mu ^2\left [ \left (\frac {R}{R_0}\right )^3 - 1 \right ] \right \}}\right )\!, \end{align}
where
$R(t) = n(t)^{-2/3}B(t)$
. In cylindrical coordinates in momentum, (2.30) reads
\begin{align} f(t,p_{\perp },p_{\parallel }) = \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + \frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^{2}(B/B_0)^{-2}}}\right )\!. \end{align}
We see that, initially, when
$n=n_0$
,
$B=B_0$
, the Maxwell–Jüttner distribution is recovered. At any later time, as long that
$n(t)$
or
$B(t)$
evolve (either growing or decreasing), a finite pressure anisotropy will develop, as we will see below. Therefore, (2.30), (2.31) constitute a natural extension of the Maxwell–Jüttner distribution to a relativistic anisotropic distribution, analogous to the bi-Maxwellian for the non-relativistic Maxwell–Boltzmann distribution.
We note that, in this case, (2.31) does not provide a natural way of obtaining expressions for
$T_{\perp }$
and
$T_{\parallel }$
directly from its functional form, as the normalised temperature
$\theta$
remains outside the square root.
By using the relativistic moment equations for the perpendicular and parallel components of the pressure tensor in a well-magnetised plasma, (cf. (2.16)), we can obtain evolution equations for
$P_{\perp }$
and
$P_{\parallel }$
in the relativistic case. For the perpendicular pressure
\begin{align} P_{\perp } &= c\int \frac {p^2(1 - \mu ^2)}{2} f \frac {d^3p}{p^0} ,\end{align}
\begin{align} P_{\perp } &= \frac {n_0c}{2 m^4c^4 \theta K_2(1/\theta )}\int \frac {1}{2}\frac {p^4(1 - \mu ^2)}{\sqrt {1 + p^2/m^2c^2}} \exp \!\left ( -\frac {1}{\theta }\sqrt {1 + A_{\mu }\frac {p^2}{m^2c^2} } \right ) {\rm d}p {\rm d}\mu , \end{align}
and for the parallel pressure
\begin{align} P_{\parallel } &= c\int p^2\mu ^2 f \frac {d^3p}{p^0} ,\end{align}
\begin{align} P_{\parallel } &= \frac {n_0c}{2 m^4c^4 \theta K_2(1/\theta )}\int \frac {p^4\mu ^2}{\sqrt {1 + p^2/m^2c^2}} \exp \!\left ( -\frac {1}{\theta }\sqrt {1 + A_{\mu }\frac {p^2}{m^2c^2} } \right ) {\rm d}p {\rm d}\mu , \end{align}
where, for simplicity, we defined
$A_{\mu } \equiv (n/n_0)^{-2/3}(R/R_0)^{-1}\{ 1 + \mu ^2[ (R/R_0)^3 - 1 ] \}$
. Note that, in the integrand of both (2.33) and (2.35), the asymmetry between the square root expressions inside the exponential term (containing the
$A_{\mu }$
term) and outside of it in the denominator makes the integral harder to solve analytically. Nevertheless, (2.33) and (2.35) provide a valid extension of the double-adiabatic equations to the relativistic regime for any temperature
$\theta$
. By numerically integrating (2.33) and (2.35), we confirm our results comparing them with the evolution of
$P_{\perp }$
and
$P_{\parallel }$
in PIC simulations of shearing and expanding boxes. Similarly, we confirm that (2.33) and (2.35) also work in the ultrarelativistic regime
$\theta \gg 1$
, validating the analytical expressions we obtain in that case (cf. (2.18), (2.21)).
3. Simulation set-up
To validate our results, we use the relativistic, 2.5-dimensional, fully kinetic PIC code TRISTAN-MP (Buneman Reference Buneman1993; Spitkovsky Reference Spitkovsky2005) to simulate a collisionless plasma made of singly charged ions and electrons subject to two types of external driving: shearing and compression. In both types of driving, ions and electrons in our simulations are initialised with Maxwell–Jüttner distributions with equal initial temperatures
$T_i^{\text{init}} = T_e^{\text{init}}$
and
$k_BT_i^{\text{init}}/m_i c^2 = 0.2$
and
$30$
. The physical parameters of our simulations are the initial temperature of ions and electrons, the initial plasma beta
$\beta _i^{\text{init}}$
, the mass ratio between ions and electrons
$m_i/m_e$
and the scale separation parameter, defined as the ratio between the non-relativistic cyclotron frequency of the ions and the shearing frequency
$\varOmega _i^{\text{init}}/s$
for the shearing simulations, and the ratio between the non-relativistic cyclotron frequency of the ions and the compressing frequency
$\varOmega _i^{\text{init}}/q$
for the compressing simulations. We chose
$\varOmega _i^{\text{init}}/s$
and
$\varOmega _i^{\text{init}}/q$
large enough to ensure adiabatic invariants are well conserved in the simulations. The numerical parameters of our simulations are the number of macroparticles per cell
$N_{\text{ppc}}$
, the plasma skin depth in terms of grid point spacing
$c/(\omega _{p,e}^2 + \omega _{p,i}^2)^{1/2}/\Delta x$
, where
$\omega _{p,e}$
and
$\omega _{p,i}$
are the electron and ion plasma frequencies, respectively, and the domain size in terms of the initial ion Larmor radius
$L/R_{L,i}$
, where
$R_{L,i} = v_{th,i}/\varOmega _i^{\text{init}}$
and
$v_{th,i}^2 = k_BT^{\text{init}}_i/m_i$
. The physical and numerical parameters of the shearing and compressing simulations presented here are listed in tables 1 and 2.
The physical and numerical parameters of the shearing simulations. The physical parameters are the initial ion plasma beta
$\beta _i^{\text{init}} = 8\pi n_ik_BT_{i}^{\text{init}}/B^2$
, the ion to electron mass ratio
$m_i/m_e$
, the initial ion temperature normalised to the rest mass energy
$k_BT_i^{\text{init}}/m_i c^2$
and the scale separation parameter
$\varOmega _i/s$
. The numerical parameters are the number of macroparticles per cell
$N_{\text{ppc}} = 300$
in all the simulations, the plasma skin depth in units of grid point spacing
$d/\Delta x \equiv c/(\omega _{p,e}^2 + \omega _{p,i}^2)^{1/2}/\Delta x$
and the domain size in units of initial ion Larmor radius
$L/R_{L,i}$
.
Table 1. Long description
The table has 8 columns and 7 rows describing physical and numerical parameters of the numerical simulations presented in the work. They are described in terms of plasma parameters and specific numerical variables.
The physical and numerical parameters of the compressing simulations. The physical parameters are the initial ion plasma beta
$\beta _i^{\text{init}} = 8\pi n_ik_BT_{i}^{\text{init}}/B^2$
, the ion to electron mass ratio
$m_i/m_e$
, the initial ion temperature normalised to the rest mass energy
$k_BT_i^{\text{init}}/m_i c^2$
and the scale separation parameter
$\varOmega _i/q$
. The numerical parameters are the number of macroparticles per cell
$N_{\text{ppc}} = 300$
in all the simulations, the plasma skin depth in units of grid point spacing
$d/\Delta x \equiv c/(\omega _{p,e}^2 + \omega _{p,i}^2)^{1/2}/\Delta x$
and the domain size in units of initial ion Larmor radius
$L/R_{L,i}$
.
Table 2. Long description
The table has 8 columns and 7 rows describing physical and numerical parameters of the numerical simulations presented in the work. They are described in terms of plasma parameters and specific numerical variables.
The shearing and compression motions will generate
$\dot {B}/B \neq 0$
and
$\dot {n}/n \neq 0$
, therefore creating a pressure anisotropy in the plasma, and
$P_{\perp ,\alpha }$
,
$P_{\parallel ,\alpha }$
(
$\alpha =i,e$
) will evolve independently. This adiabatic evolution will, in general, continue until the pressure anisotropy surpasses some instability threshold, usually scaling with
$1/\beta _{\alpha }$
to some power, at which point the pressure anisotropy will act as a free energy source to make the instability grow.
The specific instability excited will depend on the physical parameters of the plasma in each simulation. Some examples of pressure anisotropy-driven instabilities are ion mirror (Chandrasekhar, Kaufman & Watson Reference Chandrasekhar, Kaufman and Watson1958; Rudakov & Sagdeev Reference Rudakov and Sagdeev1959; Vedenov & Sagdeev Reference Vedenov and Sagdeev1959; Barnes Reference Barnes1966; Hasegawa Reference Hasegawa1969; Southwood & Kivelson Reference Southwood and Kivelson1993; Kivelson & Southwood Reference Kivelson and Southwood1996; Pokhotelov et al. Reference Pokhotelov, Treumann, Sagdeev, Balikhin, Onishchenko, Pavlenko and Sandberg2002, Reference Pokhotelov, Sagdeev, Balikhin and Treumann2004), ion firehose (Chandrasekhar et al. Reference Chandrasekhar, Kaufman and Watson1958; Parker Reference Parker1958; Vedenov & Sagdeev Reference Vedenov and Sagdeev1961; Yoon, Wu & de Assis Reference Yoon, Wu and de Assis1993; Gary et al. Reference Gary, Li, O’Rourke and Winske1998; Hellinger & Matsumoto Reference Hellinger and Matsumoto2000; Bott, Cowley & Schekochihin Reference Bott, Cowley and Schekochihin2024), ion cyclotron (Sagdeev & Shafranov Reference Sagdeev and Shafranov1960; Gary Reference Gary1992; Gary, Fuselier & Anderson Reference Gary, Fuselier and Anderson1993; López et al. Reference López, Moya, Navarro, Araneda, Muñoz, Viñas and Valdivia2016), electron mirror (Hellinger & Štverák Reference Hellinger and Štverák2018), electron firehose (Li & Habbal Reference Li and Habbal2000; Gary & Nishimura Reference Gary and Nishimura2003; López et al. Reference López, Micera, Lazar, Poedts, Lapenta, Zhukov, Boella and Shaaban2022) and electron whistler instabilities (Kennel & Petschek Reference Kennel and Petschek1966; Gary & Wang Reference Gary and Wang1996), among others (see e.g. Bott et al. Reference Bott, Cowley and Schekochihin2024 for a thorough review). Many of these instabilities have been studied in great detail using PIC simulations (e.g. Gary, Liu & Winske Reference Gary, Liu and Winske2011; Kunz, Schekochihin & Stone Reference Kunz, Schekochihin and Stone2014; Sironi & Narayan Reference Sironi and Narayan2015; Riquelme et al. Reference Riquelme, Quataert and Verscharen2015, Reference Riquelme, Quataert and Verscharen2016, Reference Riquelme, Quataert and Verscharen2018; Innocenti et al. Reference Innocenti, Tenerani, Boella and Velli2019; Bott et al. Reference Bott, Arzamasskiy, Kunz, Quataert and Squire2021; Zhdankin et al. Reference Zhdankin, Kunz and Uzdensky2023; Ley et al. Reference Ley, Zweibel, Miller and Riquelme2024).
In this work, we are interested in describing the adiabatic evolution of
$P_{\perp ,j}$
and
$P_{\parallel ,j}$
, before any instability is excited. Therefore, we will avoid capturing the excitation of any instability in the simulations by either stopping the simulation before excitation (for the shearing case), or disabling the fluctuating field evolution by not depositing particle current (for the compressing case). We show that our results are independent of the mass ratio
$m_i/m_e$
, and that they are applicable in the three initial plasma beta cases we tried:
$\beta _{i}^{\text{init}}=0.05, 0.5$
and
$5$
.
3.1. Shearing simulation set-up
For the shearing simulations, a periodic velocity shear
$\boldsymbol{v} = -sx\hat {y}$
is imposed in the domain, where
$s$
is the shearing frequency and
$x$
the distance along the
$\hat {x}$
axis. The magnetic field initially points along the
$x$
-axis,
$\boldsymbol{B} = B_0\hat {x}$
. The PIC system of equations is solved using shearing coordinates (Riquelme et al. Reference Riquelme, Quataert, Sharma and Spitkovsky2012). By flux conservation, the action of the shear continuously amplifies the magnetic field strength such that the magnitude of
$\boldsymbol{B}$
evolves as
$B(t) = B_0\sqrt {1 + s^2t^2}$
, where
$B_0$
is the initial field strength, whereas the density remains constant,
$n=n_0$
, where
$n_0$
is the initial density. Therefore, we have a contribution
$\dot {B}/B \neq 0$
and
$\dot {n}/n = 0$
. In this case, the magnetic field amplification will continuously drive an ion and electron pressure anisotropy,
$P_{\perp ,i} \gt P_{\parallel ,i}$
and
$P_{\perp ,e} \gt P_{\parallel ,e}$
, respectively, and we will see that the evolution of the perpendicular and parallel pressures will correctly be described by the relations presented in §§ 2.2 and 2.3, for mildly relativistic and ultrarelativistic initial temperatures.
3.2. Compressing simulation set-up
For the compressing simulations, a global compression motion is imposed in the domain, and the PIC system of equations is solved in comoving coordinates (Sironi & Narayan Reference Sironi and Narayan2015). To achieve this, a transformation from the laboratory frame
$(t_{\text{lab}},\boldsymbol{x}_{\text{lab}})$
to the comoving coordinate
$(t',\boldsymbol{x}')$
is performed via
$\boldsymbol{x}_{\text{lab}} = \boldsymbol{L}\boldsymbol{x}'$
, where
\begin{align} \boldsymbol{L} = \begin{pmatrix} a_x(t)\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; \\ 0\;\;\;\; & a_y(t)\;\;\;\; & 0 \\ 0\;\;\;\; & 0\;\;\;\; & a_z(t) \end{pmatrix}, \end{align}
and the differential transformation law reads
\begin{align} {\rm d}\boldsymbol{x}_{\text{lab}} = \boldsymbol{L}{\rm d}\boldsymbol{x}' + \dot {\boldsymbol{L}}\boldsymbol{x}'{\rm d}t'. \end{align}
In our two-dimensional domain, the magnetic field is aligned in the
$x$
direction, and the compression is performed along both
$y$
and
$z$
axes, i.e. along the perpendicular directions with respect to the direction of the magnetic field
$\boldsymbol{B}$
. For this, we choose the following scale factors:
\begin{align} a_x(t) &= 1, \end{align}
\begin{align} a_y(t) &= a_z(t) = \frac {1}{1 + qt}, \end{align}
where
$q$
is the compression rate, which can be controlled for each simulation. By flux freezing, both density and magnetic field strength get amplified, evolving as
\begin{align} n(t) &= n_0(1 + qt)^2 ,\end{align}
\begin{align} B(t) &= B_0(1 + qt)^2, \end{align}
where
$n_0$
is the initial density, and
$B_0$
is the initial field strength. Therefore, in this case we have
$\dot {n}/n \neq 0$
and
$\dot {B}/B \neq 0$
. The magnetic field and density amplification will also drive a continuous ion and electron pressure anisotropy,
$P_{\perp ,i} \gt P_{\parallel ,i}$
and
$P_{\perp ,e} \gt P_{\parallel ,e}$
, respectively.
It is important to note that, in the compressing simulations the parallel component of the momentum of both ions and electrons remains unchanged during the compression,
$p_{\parallel } = p_{\parallel ,0}$
. Therefore,
${\rm d}p_{\parallel }/{\rm d}t = 0$
in the drift kinetic equation (2.1). We can solve this special case following the same procedure as in § 2, and the evolution of
$P_{\perp }$
and
$P_{\parallel }$
for this case will be presented in § 4.2 along with the simulation results, also showing good agreement.
4. Simulation results
4.1. Shearing simulations
Shearing simulations generate a growing magnetic field
$B(t)/B_0 = \sqrt {1 + s^2t^2}$
and a constant density,
$n/n_0=1$
. We can then replace these relations in the relativistic double-adiabatic (2.18), (2.21) or (2.33), (2.35) depending on the initial temperature, in order to compare with simulations.
4.1.1. Initial
$k_BT_i^{\text{init}}/m_ic^2= 0.2$
The evolution of the ion perpendicular pressure
$P_{\perp ,i}$
and parallel pressure
$P_{\parallel ,i}$
as well as the evolution of the electron perpendicular and parallel pressure
$P_{\perp ,e}$
,
$P_{\parallel ,e}$
are shown in figure 1 for run Shb0.5m8d0.2wcis800.
(
$a$
) The evolution of the ion perpendicular pressure for run Shb0.5m8d0.2wcis800 (open blue circles), for initial
$k_BT_i^{\text{init}}/m_ic^2=0.2$
. The solid blue line shows the evolution of
$P_{\perp }$
according to (2.33), numerically integrated. (
$b$
) The evolution of the ion parallel pressure (open orange circles). The solid orange line shows the evolution of
$P_{\parallel }$
according to (2.35), numerically integrated. (
$c$
) The evolution of the electron perpendicular pressure (open blue circles), for initial
$k_BT_e^{\text{init}}/m_ec^2=1.6$
(given a mass ratio
$m_i/m_e = 8$
). The solid blue line shows the evolution of
$P_{\perp }$
according to (2.33), numerically integrated. (
$d$
) The evolution of the electron parallel pressure (open orange circles). The solid orange line shows the evolution of
$P_{\parallel }$
according to (2.35), numerically integrated. In panels (
$a$
) and (
$c$
), the non-relativistic CGL evolution for
$P_{\perp }$
is shown in dashed grey line. In panels (
$b$
) and (
$d$
), the non-relativistic CGL evolution for
$P_{\parallel }$
is shown in dashed grey line.
Figure 1. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the solid lines fully coincide with the open circles lines, and the dashed line agree only at the beginning, and diverges at later times. Read caption for more details.
For this mildly relativistic temperature, we can already see that
$P_{\perp ,i}$
(open blue circles in panel
$a$
) departs from the non-relativistic CGL prediction for a pure shearing,
$B/B_0$
(dashed grey line in panel
$a$
, cf. (2.11)). In the parallel pressure (open orange circles in panel
$b$
), we also see an analogous departure from the CGL prediction
$(B/B_0)^{-2}$
, although it is less pronounced than in the perpendicular case.
Our relativistic extension for the evolution of
$P_{\perp }$
is shown in the solid blue line in panel
$a$
of figure 1. This line is obtained by numerically integrating (2.33) for
$\theta =0.2$
. We can see that it consistently follows the simulation results for
$P_{\perp ,i}$
. Similarly, the relativistic extension for the evolution of
$P_{\parallel }$
, obtained by numerically integrating (2.35) for
$\theta =0.2$
and shown in the solid orange line in panel
$b$
, also shows very good agreement with the evolution of
$P_{\parallel ,i}$
. We can see that our derivations outlined in § 2.3 correctly describe the behaviour of
$P_{\perp }$
and
$P_{\parallel }$
for this mildly relativistic temperature.
We see an analogous behaviour for electrons. Given that
$m_i/m_e=8$
, we initially have
$k_BT_e^{\text{init}}/m_ec^2 = k_BT_i^{\text{init}}/m_ic^2 (m_i/m_e) = 1.6$
. By numerically integrating (2.33), (2.35) for
$\theta =1.6$
, we also obtain an evolution of
$P_{\perp }$
and
$P_{\parallel }$
that is consistent with our PIC results. Indeed, we can see the evolution of
$P_{\perp ,e}$
(open blue circles, panel
$c$
) and
$P_{\parallel ,e}$
(open orange circles, panel
$d$
) for run Shb0.5m8d0.2wcis800. Equations (2.33) (solid blue line in panel
$c$
) and (2.35) (solid orange line in panel
$d$
) correctly describe the evolution of
$P_{\perp ,e}$
and
$P_{\parallel ,e}$
, respectively. In this case, the departure from the non-relativistic CGL evolution is more pronounced in both perpendicular and parallel pressures, as the initial temperature of electrons is higher by a factor
$m_i/m_e=8$
.
4.1.2. Initial
$k_BT_i^{\text{init}}/m_ic^2 = 30$
The evolution of the ion perpendicular pressure
$P_{\perp ,i}$
and parallel pressure
$P_{\parallel ,i}$
, as well as the evolution of the electron perpendicular and parallel pressures
$P_{\perp ,e}$
,
$P_{\parallel ,e}$
are shown in figure 2 for run Shb0.5m1836d30wcis3200. For initial
$k_BT_i^{\text{init}}/m_ic^2 = 30$
, we are already in the ultrarelativistic regime, so we can now test our analytic derivation of equations (2.18) and (2.21).
For the ion perpendicular pressure
$P_{\perp ,i}$
(panel
$a$
, open black circles), we can see that it departs from the non-relativistic CGL prediction
$B/B_0$
(dashed grey line) much further than in the lower initial temperature case. In contrast, we observe that our ultrarelativistic extension for the evolution of
$P_{\perp }$
, equation (2.18) (dashed blue line in panel
$a$
), consistently follows the evolution of
$P_{\perp ,i}$
. Our general equation (2.33) (numerically integrated for
$\theta =30$
, solid orange line) is also shown, and it also consistently follows the evolution of
$P_{\perp ,i}$
. Therefore, our analytical expression for
$P_{\perp }$
in the ultrarelativistic case is consistent with our general formula, equation (2.33).
(
$a$
) The evolution of the ion perpendicular pressure for run Shb0.5m1836d30wcis3200 (open black circles), for initial
$k_BT_i^{\text{init}}/m_ic^2=30$
. The dashed blue line shows the evolution of
$P_{\perp }$
according to equation (2.18). The solid orange line shows the evolution of
$P_{\perp }$
according to equation (2.33), numerically integrated. (
$b$
) The evolution of the ion parallel pressure (open black circles). The dashed magenta line shows the evolution of
$P_{\parallel }$
according to equation (2.21). The solid green line shows the evolution of
$P_{\parallel }$
according to equation (2.35), numerically integrated. (
$c$
) The evolution of the electron perpendicular pressure (open black circles), for initial
$k_BT_e^{\text{init}}/m_ec^2=55\,080$
. The dashed blue line shows the evolution of
$P_{\perp }$
according to (2.18). The solid orange line shows the evolution of
$P_{\perp }$
according to (2.33), numerically integrated. (
$d$
) The evolution of the electron parallel pressure (open black circles). The dashed magenta line shows the evolution of
$P_{\parallel }$
according to (2.21). The solid green line shows the evolution of
$P_{\parallel }$
according to (2.35), numerically integrated. In panels
$(a)$
and
$(c)$
, the non-relativistic CGL evolution for
$P_{\perp }$
is shown in dashed grey line. In panels
$(b)$
and
$(d)$
, the non-relativistic CGL evolution for
$P_{\parallel }$
is shown in dashed grey line.
Figure 2. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the colored solid lines fully coincide with the open circles lines, and the gray dashed line agree only at the beginning, and diverges at later times. Read caption for more details.
The evolution of the ion parallel pressure
$P_{\parallel ,i}$
(open black circles, right panel of figure 2) follows a similar trend. The departure from the non-relativistic CGL prediction
$(B/B_0)^{-2}$
(dashed grey line) is also evident, and our ultrarelativistic extension for
$P_{\parallel }$
, equation (2.21) (solid magenta line in figure 2), is also consistent with the simulation results. We show the general formula (2.35) (numerically integrated for
$\theta =30$
, solid green line), which also explains very well the evolution of
$P_{\parallel ,i}$
. Therefore, our ultrarelativistic analytical result is also consistent with the general formula, equation (2.35).
Similar to the previous case, the electron pressure in run Shb0.5m1836d30wcis3200 has an analogous evolution, now with a much higher initial temperature, given the realistic mass ratio
$m_i/m_e=1836$
. Nevertheless, the ultrarelativistic extension derived in § 2.2 is also valid, and correctly describes the evolution of the electron pressure components.
Indeed, the evolution of the electron perpendicular pressure
$P_{\perp ,e}$
and parallel pressure
$P_{\parallel ,e}$
for run Shb0.5m1836d30wcis3200 is shown in figure 2. Similar to the ion case, we can see that neither
$P_{\perp ,e}$
(panel
$c$
, open black circles) nor
$P_{\parallel ,e}$
(panel
$d$
, open black circles) follow the respective CGL prediction (dashed grey lines in panels
$c$
and
$d$
of figure 2), whereas our ultrarelativistic extensions, equations (2.18) (dashed blue line in panel
$c$
) and (2.21) (dashed magenta line in panel
$d$
) correctly describes their evolution. We also include the general equation (2.33) (solid orange line, panel
$c$
) and equation (2.35) (solid green line, panel
$d$
) numerically integrated for
$\theta =30$
, showing good agreement as well.
Our ultrarelativistic extensions of the double-adiabatic evolutions of
$P_{\perp }$
, equation (2.18) and
$P_{\parallel }$
, equation (2.21) are demonstrated to agree with PIC simulations with a shear external driving. Additionally, they do not depend on the initial temperature, provided that it is sufficiently high,Footnote
3
analogous to the non-relativistic, CGL case, which is also independent of temperature (as long as temperatures are sufficiently low).
4.2. Compressing simulations
As mentioned in § 3, the compressing simulations leave the parallel momentum unchanged,
$p_{\parallel }=p_{\parallel ,0}$
. The constancy of parallel momentum implies that
${\rm d}p_{\parallel }/{\rm d}t = 0$
in the drift kinetic equation (2.1)
\begin{align} \frac {\partial f}{\partial t} + \frac {{\rm d}p_{\perp }}{{\rm d}t}\frac {\partial f}{\partial p_{\perp }} = \frac {\partial f}{\partial t} + \frac {\dot {B}}{2B}p_{\perp }\frac {\partial f}{\partial p_{\perp }} = 0. \end{align}
This just constitutes a special case of the general derivation presented in § 2. We note that in this case the drift kinetic equation can also readily be solved by the method of characteristics to obtain
\begin{align} f(t,p_{\perp },p_{\parallel }) = f_0\left (\frac {p_{\perp }}{(B/B_0)^{1/2}}, p_{\parallel }\right )\!. \end{align}
For an ultrarelativistic Maxwell–Jüttner distribution as initial condition, our solution (4.2) becomes
\begin{align} f(t,p_{\perp },p_{\parallel }) = \frac {n_0}{8\pi p_T^3}\exp \!\left (-\frac {1}{p_T}\sqrt {\frac {p_{\perp }^2}{B/B_0} + p_{\parallel }^2}\right )\!, \end{align}
and taking the perpendicular and parallel pressure moments (cf. eqns. (2.16))
\begin{align} & P_{\perp } = \frac {3}{4}P_0\left (\frac {B}{B_0}\right )^2\left [ \frac {(B/B_0)-2}{((B/B_0)-1)^{3/2}}\arctan \left (\sqrt {(B/B_0)-1}\right ) + \frac {1}{(B/B_0)-1} \right ]\\& \nonumber \end{align}
for
$(B/B_0)\gt 1$
, and
\begin{align} P_{\perp } = \frac {3}{4}P_0\left (\frac {B}{B_0}\right )^2\left [ \frac {2 - (B/B_0)}{(1-(B/B_0))^{3/2}}\text{arctanh}\left ( \sqrt {1 - (B/B_0)} \right ) - \frac {1}{1 - (B/B_0)} \right ] \end{align}
for
$(B/B_0)\lt 1$
. Similarly for the parallel pressure
\begin{align} P_{\parallel } = \frac {3}{2}P_0\left [ \frac {(B/B_0)^2}{((B/B_0) - 1)^{3/2}}\arctan \left ( \sqrt {(B/B_0)-1} \right ) - \frac {(B/B_0)}{(B/B_0)-1} \right ] \end{align}
for
$(B/B_0)\gt 1$
, and
\begin{align} P_{\parallel } = \frac {3}{2}P_0\left [ \frac {(B/B_0)}{1 - (B/B_0)} - \frac {(B/B_0)^2}{(1 - (B/B_0))^{3/2}}\text{arctanh}\left (\sqrt {1 - (B/B_0)}\right ) \right ] \end{align}
for
$(B/B_0)\lt 1$
.
For an initial Maxwell–Jüttner distribution, our solution (4.2) becomes
\begin{align} f(t,p_{\perp },p_{\parallel }) = \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\exp \!\left ( -\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + p_{\parallel }^2/m^2c^2} \right )\!, \end{align}
and the perpendicular and parallel pressures can also be obtained by numerically integrating the moments
\begin{align} P_{\perp } &= \frac {n_0c}{2m^4c^4\theta K_2(1/\theta )} \int \frac {1}{2}\frac {p_{\perp }^2}{\sqrt {1 + p_{\perp }^2/m^2c^2 + p_{\parallel }^2/m^2c^2}}\nonumber \\[4pt] &\quad \exp \!\left ( -\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + p_{\parallel }^2/m^2c^2} \right )p_{\perp }{\rm d}p_{\perp }{\rm d}p_{\parallel } , \end{align}
and
\begin{align} P_{\parallel } &= \frac {n_0c}{2m^4c^4\theta K_2(1/\theta )} \int \frac {p_{\parallel }^2}{\sqrt {1 + p_{\perp }^2/m^2c^2 + p_{\parallel }^2/m^2c^2}}\nonumber \\[4pt] &\quad \exp \!\left ( -\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + p_{\parallel }^2/m^2c^2} \right )p_{\perp }{\rm d}p_{\perp }{\rm d}p_{\parallel }. \end{align}
We can now use equations (4.4), (4.6), (4.9) and (4.10) to compare with our compressing box simulations. In all the results shown below, we pass from pressures to temperatures using
$T_{\perp }=P_{\perp }/(n k_B)$
and
$T_{\parallel } = P_{\parallel }/(n k_B)$
, where
$n(t)=n_0(1 + qt)^2$
, as in the relativistic regime it also holds that
$P=nk_BT$
(e.g. Rezzolla & Zanotti Reference Rezzolla and Zanotti2013).
4.2.1. Initial
$k_BT_i^{\text{init}}/m_i c^2 = 0.2$
(
$a$
) The evolution of the ion perpendicular temperature for run Compb0.5m8d0.2wcis800 (open blue circles), for initial
$k_BT_i^{\text{init}}/m_ic^2=0.2$
. The solid blue line shows the evolution of
$T_{\perp }$
according to (4.9), numerically integrated. (
$b$
) The evolution of the ion parallel temperature (open orange circles). The solid orange line shows the evolution of
$T_{\parallel }$
according to (4.10), numerically integrated. (
$c$
) The evolution of the electron perpendicular temperature (open blue circles), for initial
$k_BT_e^{\text{init}}/m_ec^2=1.6$
. The solid blue line shows the evolution of
$T_{\perp }$
according to (4.9), numerically integrated. Panel
$(d)$
the evolution of the electron parallel temperature (open orange circles). The solid orange line shows the evolution of
$T_{\parallel }$
according to (4.10), numerically integrated. In panels
$(a)$
and
$(c)$
the non-relativistic CGL evolution for
$P_{\perp }$
is shown in dashed grey line for a compressing motion. In panels
$(b)$
and
$(d)$
, the non-relativistic CGL evolution for
$P_{\parallel }$
is shown in dashed grey line for a compressing motion.
Figure 3. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the solid lines fully coincide with the open circles lines, and the dashed line agree only at the beginning, and diverges at later times. Read caption for more details.
The evolutions of the ion perpendicular temperature
$T_{\perp ,i}$
and parallel temperature
$T_{\parallel ,i}$
, and the electron perpendicular and parallel temperatures
$T_{\perp ,e}$
,
$T_{\parallel ,e}$
are shown in figure 3 for run Compb0.5m8d0.2wcis800.
Similar to the shearing case, we can see that both
$T_{\perp ,i}$
(open blue circles in panel
$a$
) and
$T_{\parallel ,i}$
(open orange circles in panel
$b$
) depart from the non-relativistic CGL prediction (dashed grey lines). This is especially critical in the parallel case, as the non-relativistic CGL predicts a constant
$T_{\parallel }$
. This departure from CGL was also observed in Tran et al. (Reference Tran, Sironi, Ley, Zweibel and Riquelme2023) for the same initial ion temperature.
In contrast, we can see that our relativistic double-adiabatic equations describe very well the temperature evolution in our simulations. We numerically integrated equations (4.9) (solid blue line in panel
$a$
) and (4.10) (solid orange line in panel
$b$
) for
$\theta = 0.2$
.
The evolution of the electron perpendicular temperature
$T_{\perp ,e}$
and parallel temperature
$T_{\parallel ,e}$
for run Compb0.5m8d0.2wcis800 is shown in figure 3. In this case we also see a very similar behaviour, now for
$k_BT_e^{\text{init}}=1.6$
. Both
$T_{\perp ,e}$
(open blue circles, panel
$c$
) and
$T_{\parallel ,e}$
(open orange circles, panel
$d$
) depart from their non-relativistic CGL evolution (dashed grey lines), and are very well described by our relativistic double-adiabatic equations (4.9), (4.10).
4.2.2. Initial
$k_BT_i^{\text{init}}/m_ic^2 = 30$
The evolution of the ion perpendicular temperature
$T_{\perp ,i}$
and parallel temperature
$T_{\parallel ,i}$
, and the electron perpendicular and parallel temperatures
$T_{\perp ,e}$
,
$T_{\parallel ,e}$
are shown in figure 4 for run Compb0.5m8d0.2wcis800. Here, we can test our ultrarelativistic double-adiabatic equations (4.4) and (4.6).
(
$a$
) The evolution of the ion perpendicular temperature for run Compb0.5m1836d30wcis3200 (open blue circles), for initial
$k_BT_i^{\text{init}}/m_ic^2=30$
. The solid blue line shows the evolution of
$T_{\perp }$
according to (4.4). (
$b$
) The evolution of the ion parallel temperature (open orange circles). The solid orange line shows the evolution of
$T_{\parallel }$
according to (4.6). In panels
$(a)$
and
$(c)$
the non-relativistic CGL evolution for
$P_{\perp }$
is shown in dashed grey line for a compressing motion. In panels
$(b)$
and
$(d)$
, the non-relativistic CGL evolution for
$P_{\parallel }$
is shown in dashed grey line for a compressing motion.
Figure 4. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the solid lines fully coincide with the open circles lines, and the dashed line do not follow the rest of the lines at all. Read caption for more details.
Similar to the shearing case, we can see that, for this ultrarelativistic initial temperature, both
$T_{\perp ,i}$
(open blue circles in panel
$a$
) and
$T_{\parallel ,i}$
(open orange circles in panel
$b$
) also depart from their non-relativistic CGL prediction (dashed grey lines), especially for
$T_{\perp ,i}$
, which shows a more pronounced deviation than in the mildly relativistic case, due to the high initial temperature.
Nevertheless, we also see here that our ultrarelativistic double-adiabatic equations (4.4) (solid blue line in panel
$a$
) and (4.6) (solid orange line in panel
$b$
) consistently follow the evolution of
$T_{\perp ,i}$
and
$T_{\parallel ,i}$
from the simulation.
The evolution of the electron perpendicular temperature
$T_{\perp ,e}$
and parallel temperature
$T_{\parallel ,e}$
for run Compb0.5m8d0.2wcis800 is shown in figure 4.
For electrons in run Compb0.5m1836d30wcis3200, the behaviour is analogous to the corresponding shearing simulation. We see that the departure from non-relativistic CGL predictions (dashed grey lines) is present in both
$P_{\perp ,e}$
(open blue circles, panel
$c$
) and
$P_{\parallel ,e}$
(open orange circles, panel
$d$
). For this very high
$k_BT_e^{\text{init}}/m_ec^2$
, we can see that our ultrarelativistic double-adiabatic equations (4.4) (solid blue line, panel
$c$
), (4.6) (solid orange line, panel
$d$
) also work very well in describing the evolution of
$P_{\perp ,e}$
and
$P_{\parallel ,e}$
in the simulations.
4.3. Distribution function comparison
From both shearing and compressing simulations, we can obtain the ion and electron distribution functions at any time by constructing the histograms of the particle momenta
$p_{\perp ,\alpha }$
,
$p_{\parallel ,\alpha }$
(
$\alpha =i,e$
). We can then compare these empirical distributions with our solution (2.31) of the drift kinetic equation (2.2), to assess how well it describes the evolution of the distribution in the simulations, from an initial, isotropic Maxwell–Jüttner distribution into its anisotropic form.
From our time-dependent distribution function, equation (2.31),
\begin{align}& f(t,p_{\perp },p_{\parallel }) = \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + \frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^{2}(B/B_0)^{-2}}}\right )\!,\nonumber \\& \end{align}
we can obtain the marginal distributions in
$p_{\perp }$
and
$p_{\parallel }$
by integrating
\begin{align} f_{\text{aMJ}}(t,p_{\perp }) &= \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )} \int _{-\infty }^{\infty }\nonumber \\[5pt] &\quad \exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + \frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^{2}(B/B_0)^{-2}}}\right ) 2\pi p_{\perp } {\rm d}p_{\parallel },\\[-8pt] &\nonumber \end{align}
\begin{align} f_{\text{aMJ}}(t,p_{\perp }) &= \frac {n_0}{m^2c^2 \theta K_2(1/\theta )} \frac {(n/n_0)}{(B/B_0)}p_{\perp }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)}} K_1\left (\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)}}\right )\!,\nonumber \\&\end{align}
and
\begin{align}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! f_{\text{aMJ}}(t,p_{\parallel }) &= \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )} \int _{0}^{\infty }\nonumber \\[5pt] &\quad \exp \!\left (-\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + \frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^{2}(B/B_0)^{-2}}}\right ) 2\pi p_{\perp } {\rm d}p_{\perp } ,\end{align}
\begin{align}\,\, f_{\text{aMJ}}(t,\!p_{\parallel }) &= \frac {n_0}{2mc K_2(1/\theta )}\!\left (\!\sqrt {1\! + \!\frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^2(B/B_0)^{-2}}} + \theta \!\right )\! \exp \!\left ( -\frac {1}{\theta }\sqrt {\!1\! +\! \frac {p_{\parallel }^2/m^2c^2}{(n/n_0)^2(B/B_0)^{-2}}} \right )\!,\end{align}
where
$K_1(x)$
is the modified Bessel function of the first kind. The comparison between the histograms of ion perpendicular momentum
$p_{\perp ,i}$
with the marginal distribution
$f_{\text{aMJ}}(t,p_{\perp })$
and the ion parallel momentum
$p_{\parallel ,i}$
with
$f_{\text{aMJ}}(t,p_{\parallel })$
for our shearing simulations is shown in figure 5 for initial
$k_BT_i^{\text{init}}/m_ic^2 = 0.2$
and
$k_BT_i^{\text{init}}/m_ic^2 = 30$
. This is shown at the beginning and at the end of the simulations (i.e. from initial isotropic state at
$t\boldsymbol{\cdot }s=0$
to an anisotropic state at
$t\boldsymbol{\cdot }s= 1$
). We are showing ion results only, as electron results show a very similar behaviour.
(
$a$
) The histograms of the ion perpendicular momentum for run Shb0.5m8d0.2wcis800 at
$t\boldsymbol{\cdot }s = 0$
(light blue bars) and
$t\boldsymbol{\cdot }s = 1$
(yellow bars). The marginal distribution
$f_{\text{aMJ}}(p_{\perp },t)$
, (4.13), is shown at
$t\boldsymbol{\cdot }s = 0$
(solid blue line) and
$t\boldsymbol{\cdot }s = 1$
(solid red line), for
$\theta =0.2$
. (
$b$
) The histograms of the ion parallel momentum at
$t\boldsymbol{\cdot }s = 0$
(light blue bars) and
$t\boldsymbol{\cdot }s = 1$
(yellow bars). The marginal distribution
$f_{\text{aMJ}}(p_{\parallel },t)$
, (4.15), is shown at
$t\boldsymbol{\cdot }s = 0$
(solid blue line) and
$t\boldsymbol{\cdot }s = 1$
(solid red line), for
$\theta =0.2$
. Panels
$(c)$
and
$(d)$
show the same quantities as panels
$(a)$
and
$(b)$
, respectively, but for run Shb0.5m8d30wcis3200 (i.e.
$\theta =30$
).
Figure 5. Long description
A mosaic of four graphs describing the histogram of specific variables and a theoretical prediction for the distribution function in those variables, at two different times in a specific simulation run. In each of the four panels, the solid lines fully coincide with the envelope of the histograms for each time in the simulation that is presented. Read caption for more details.
We can see that, for all cases presented, the agreement between simulations and theory is noteworthy. For the perpendicular momentum (panels
$a$
and
$c$
in figure 5) we can see very good agreement at
$t\boldsymbol{\cdot }s = 0$
between the
$p_{\perp ,i}$
histogram (light blue bars) and our analytical marginal distribution
$f_{\text{aMJ}}(ts=0,p_{\perp })$
((4.13), solid blue line), validating the isotropic limit to the Maxwell–Jüttner distribution. More notably, at the end of the simulation we see the expected widening of the
$p_{\perp ,i}$
histogram (yellow bars) towards higher values, transitioning to an anisotropic state, and our marginal distribution
$f_{\text{aMJ}}(ts=1,p_{\perp })$
(solid red line) shows very good agreement with the simulation at this stage.
For the parallel momentum (panels
$b$
and
$d$
in figure 5) we also obtain very good correspondence at
$t\boldsymbol{\cdot }s = 0$
between
$p_{\parallel ,i}$
(light blue bars) and the analytical marginal distribution
$f_{\text{aMJ}}(ts=0,p_{\parallel })$
((4.15), solid blue line), and we validate the isotropic limit in this case as well. At
$t\boldsymbol{\cdot }s = 1$
, we see the expected narrowing of the parallel momentum distribution (yellow bars) towards smaller values, and in this case we see that our marginal distribution
$f_{\text{aMJ}}(ts=1,p_{\parallel })$
(solid red line) correctly characterises this anisotropic state as well.
(
$a$
) The histograms of the ion perpendicular momentum for run Compb0.5m8d0.2wcis800 at
$t\boldsymbol{\cdot }s = 0$
(light blue bars) and
$t\boldsymbol{\cdot }s = 1$
(yellow bars). The marginal distribution
$f_{\text{aMJ,c}}(p_{\perp },t)$
, (4.18), is shown at
$t\boldsymbol{\cdot }s = 0$
(solid blue line) and
$t\boldsymbol{\cdot }s = 1$
(solid red line), for
$\theta =0.2$
. (
$b$
) The histograms of the ion perpendicular momentum for run Compb0.5m8d30wcis3200 at
$t\boldsymbol{\cdot }s = 0$
(light blue bars) and
$t\boldsymbol{\cdot }s = 1$
(yellow bars). The marginal distribution
$f_{\text{aMJ,c}}(p_{\perp },t)$
, (4.18), is shown at
$t\boldsymbol{\cdot }s = 0$
(solid blue line) and
$t\boldsymbol{\cdot }s = 1$
(solid red line), for
$\theta =30$
.
Figure 6. Long description
A mosaic of two graphs describing the histogram of specific variables and a theoretical prediction for the distribution function in those variables, at two different times in a specific simulation run. In each of the two panels, the solid lines fully coincide with the envelope of the histograms for each time in the simulation that is presented. Read caption for more details.
For our compressing simulations, recall the solution (4.8)
\begin{align} f(t,p_{\perp },p_{\parallel }) = \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )}\exp \!\left ( -\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + p_{\parallel }^2/m^2c^2} \right )\!. \end{align}
The marginal distribution
$f_{\text{aMJ,c}}(t,p_{\perp })$
can be calculated as
\begin{align} f_{\text{aMJ,c}}(t,p_{\perp }) &= \frac {n_0}{4\pi m^3c^3 \theta K_2(1/\theta )} \int _{-\infty }^{\infty } \!\exp \!\left (\!-\frac {1}{\theta }\sqrt {1 + \frac {p_{\perp }^2/m^2c^2}{(B/B_0)} + p_{\parallel }^2/m^2c^2}\right )\! 2\pi p_{\perp } {\rm d}p_{\parallel } ,\end{align}
The comparison between the histograms of ion perpendicular momentum
$p_{\perp ,i}$
with
$f_{\text{aMJ,c}}(t,p_{\perp })$
for the compressing simulations is shown in figure 6, also for
$k_BT_i^{\text{init}}/m_ic^2 = 0.2$
and
$k_BT_i^{\text{init}}/m_ic^2 = 30$
, and at
$t\boldsymbol{\cdot }q = 0$
and
$t\boldsymbol{\cdot }q = 1$
. We only show the perpendicular momentum, as we know that the parallel momentum of particles does not evolve. Electron results are not shown but are completely analogous.
The agreement between
$f_{\text{aMJ,c}}(t,p_{\perp })$
and
$p_{\perp ,i}$
from simulations is also remarkable in this case, for both initial temperature cases. At
$t\boldsymbol{\cdot }q = 0$
, the histogram of the perpendicular momentum (light blue bars) is very well described by
$f_{\text{aMJ,c}}(tq=0,p_{\perp })$
(solid blue line), validating the isotropic limit of our solution in this case as well. At the end of the simulation, we can clearly see the expected spread of the
$p_{\perp ,i}$
distribution (yellow bars) towards higher values, and the marginal
$f_{\text{aMJ,c}}(tq=1,p_{\perp })$
(solid red line) describes this anisotropic state very accurately as well.
We have therefore demonstrated that the distribution function (2.30), (2.31), which is a solution to the drift kinetic equation (2.2), correctly describes the double-adiabatic evolution of the distribution in shearing and compressing PIC simulations at any time. Therefore, these distributions constitute a newly proposed form for a time-dependent, anisotropic Maxwell–Jüttner distribution, derived from first principles.
5. Discussion
The results reported in this work are relevant for a diversity of high-energy astrophysical phenomena involving relativistically hot, collisionless plasmas, at both low and high beta. Among these systems one can mention pulsar wind nebulae, low-luminosity accretion flows around SMBH, relativistic jets from SMBH in active galaxies like blazars and also cosmic rays. In absence of heat fluxes and where spatial gradients are much larger than the gyroradii of particles, the relativistic species present in these plasmas can be modelled by the distribution functions (2.30), (2.31) and (2.14) for ultrarelativistic temperatures. Additionally, the adiabatic evolution of the perpendicular and parallel pressures associated with these relativistic species are well described by equations (2.33), (2.35) and (2.18), (2.21). This evolution will be valid until kinetic microinstabilities are excited and start to provide an effective collisionality in the plasma via complex wave–particle interactions.
5.1. Relativistic extended MHD in the Braginskii limit
A direct application of these results is in the context of extended MHD theories for relativistic collisionless plasmas (e.g. TenBarge, Hazeltine & Mahajan Reference TenBarge, Hazeltine and Mahajan2008; Chandra et al. Reference Chandra, Gammie, Foucart and Quataert2015). The distribution function presented in this work, along with the relativistic double-adiabatic equations, can be used as a closure to the MHD fluid system while retaining pressure anisotropies in the stress–energy tensor, provided that the system has no heat fluxes and spatial gradients in
$B$
and
$n$
are much larger than the particles’ gyration radius.
Interestingly, we can extend the results in Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015) for the evolution of the pressure anisotropy
$\Delta P = P_{\perp }-P_{\parallel }$
in the ultrarelativistic case, and apply their framework to the collisional Braginskii limit (Braginskii Reference Braginskii1965). Considering the
$R\gt 1$
case only, from equations (2.18), (2.21) and
$P=(2P_{\perp } + P_{\parallel })/3$
we can write
\begin{align} \frac {1}{P}\frac {{\rm d}\Delta P}{{\rm d}t} &= \mathcal{A}\frac {\dot {B}}{B} + \mathcal{C}\frac {\dot {n}}{n}, \end{align}
where
\begin{align} \mathcal{A} &= \frac {3}{2}\varPsi (r)\left [ \frac {(r^3 + 2)(r^3 + 8)}{2(r^3 - 1)^{5/2}}\arctan \left ( \sqrt {r^3 - 1} \right ) + \frac {r^6 - 32r^3 + 4}{2r^3(r^3-1)^2} \right ], \end{align}
\begin{align} \mathcal{C} &= -3\varPsi (r)\left [ \frac {r^3(10-r^3)}{2(r^3 - 1)^{5/2}}\arctan \left(\sqrt {r^3 - 1}\right) - \frac {r^6 + 12r^3 - 4}{2r^3(r^3-1)^2} \right ], \end{align}
\begin{align} \varPsi (r) &= \frac {r^3\sqrt {r^3 - 1}}{r^3\arctan (\sqrt {r^3 - 1}) + \sqrt {r^3 - 1}}, \end{align}
and
$r \equiv R/R_0$
. Equation (5.1) is valid for any value of anisotropy (as long as the evolution is adiabatic, see below). Note that, when we go to the small anisotropy limit, corresponding to
$r \rightarrow 1$
(i.e.
$n$
and
$B$
close to their initial values), we exactly recover the result in Chandra et al. (Reference Chandra, Gammie, Foucart and Quataert2015). Indeed, as
$\varPsi (r\rightarrow 1)\rightarrow 1/2$
,
$\mathcal{A} \rightarrow 12/5$
and
$\mathcal{C}\rightarrow 8/5$
, we get
\begin{align} \frac {1}{P}\frac {{\rm d}\Delta P}{{\rm d}t} = \frac {12}{5}\frac {\dot {B}}{B} - \frac {8}{5}\frac {\dot {n}}{n} = \frac {4}{5}\frac {d}{{\rm d}t}\ln \left (\frac {B^3}{n^2}\right ) + \mathcal{O}\left ([\Delta P/P]^2\right )\!. \end{align}
In the Braginskii limit, the adiabatic evolution of the pressure anisotropy is balanced with the scattering process present in the plasma. Denoting the scattering rate as
$\nu$
in the collisional limit, we can write, in equilibrium,
\begin{align} \Delta P = \frac {P}{\nu }\left [ \mathcal{A}\frac {\dot {B}}{B} - \mathcal{C}\frac {\dot {n}}{n} \right ]. \end{align}
Equation (5.6) is then valid for any values of
$B(t)$
and
$n(t)$
, in the ultrarelativistic regime (
$k_BT/mc^2 \gtrsim 1$
, in practice).
5.2. Small-amplitude limit
The analysis up to now is valid for arbitrarily large changes in
$B$
and
$n$
, provided the changes are slow and the plasma is collisionless, due to the isotropising effect of collisions. However, even when binary collisions are negligible, particles can be scattered by interactions with waves, and a case of particular interest is scattering by waves generated by the particles themselves, through kinetic scale anisotropy instabilities. When the threshold for these instabilities is low, the evolution up to the triggering of instability can be calculated from perturbation theory. For example, cosmic rays trigger a gyroresonant instability of Alfvén waves when the fractional pressure anisotropy is of order
$v_A/c$
(e.g. Zweibel Reference Zweibel2020).
Suppose we replace
$R/R_0$
and
$n/n_0$
by
$1+\delta R/R$
and
$1+\delta n/(3n)$
, respectively, where
$\vert \delta R/R\vert$
and
$\vert \delta n/n\vert$
are
$\ll 1$
. And assume that
$f_0$
is isotropic. Then (2.6) can be replaced by a Taylor expansion, which to first order in the perturbation amplitude is
\begin{align} f(p,\mu )=f_0(p)+\left ( \frac {\delta R}{R}P_2(\mu )-\frac {\delta n}{3n}P_0(\mu ) \right )p\frac {{\rm d}f_0}{{\rm d}p}, \end{align}
where
$P_2$
and
$P_0$
are Legendre functions of order 2 and 0, respectively.
Equation (5.7) is in convenient form for calculation of the pressure anisotropy
\begin{align} \Delta P\equiv P_{\perp } - P_{\parallel }=-\int p^2{\rm d}p {\rm d}\mu pvP_2(\mu )f(p,\mu ). \end{align}
\begin{align} \Delta P =-\frac {2}{5}\frac {\delta R}{R}\int p^4v\frac {{\rm d}f_0}{{\rm d}p}. \end{align}
For the important example of cosmic rays, we take the ultrarelativistic limit
$v\rightarrow c$
. Then, integrating (5.9) by parts and using the standard definition of pressure, we find
\begin{align} \Delta P =\frac {12}{5}\frac {\delta R}{R}P. \end{align}
6. Conclusions
In this work, we have extended the double-adiabatic equations for a homogeneous, collisionless plasma to the relativistic regime, and presented a novel anisotropic version of the Maxwell–Jüttner distribution. We analytically solved the drift kinetic equation in the homogeneous, collisionless limit, and obtained a general, time-dependent solution for the distribution function in terms of density and magnetic field variations.
By considering a relativistic Maxwell–Jüttner distribution as initial condition, we obtained new, relativistic expressions for the evolution of both perpendicular and parallel pressures, in terms of the integral moments that can be solved numerically for any temperature
$\theta = k_BT/mc^2$
. When considering an ultrarelativistic Maxwell–Jüttner distribution as initial condition, the integral moments are tractable, and we obtained analytical evolution equations for the perpendicular and parallel pressures. In practice, these analytical expressions prove to be applicable for
$k_BT/mc^2 \geq 1$
. We confirmed that our relativistic double-adiabatic equations agree with the general state equations for a relativistic anisotropic plasma from Gedalin (Reference Gedalin1991), and we validated them using PIC simulations of shearing and compressing boxes for mildly relativistic and ultrarelativistic initial temperatures, obtaining remarkable agreement. From a general perspective, we conclude that the specific form of the double-adiabatic equations then depends on the initial distribution function that is used.
Additionally, the time-dependent solution (2.31) obtained for an initial Maxwell–Jüttner distribution constitutes a new anisotropic extension for a relativistic Maxwellian, explicitly dependent on the density and magnetic field variations. For a specific value of the pressure anisotropy, the corresponding double-adiabatic equations can be used to map the density and magnetic field variations to the specific level of anisotropy. We confirmed that this distribution correctly describes the momentum distributions of both shearing and compressing PIC simulations, in both initial isotropic and final anisotropic states. As PIC simulations utilise a discretised version of the method of characteristics in their algorithms, we argue that our distribution (2.31) is the analytical form of the distribution function that our shearing and compressing box simulations evolve (when the initial condition is a Maxwell–Jüttner distribution), until an instability is excited. This analytical anisotropic Maxwell–Jüttner form could also be used to qualitatively assess the distortions to the numerical distribution function as a consequence of wave–particle interactions by microinstabilities in driven PIC simulations.
These new findings open several paths for future work. The relativistic double-adiabatic equations can be used as a new closure for extended models of relativistic MHD equations (TenBarge et al. Reference TenBarge, Hazeltine and Mahajan2008; Chandra et al. Reference Chandra, Gammie, Foucart and Quataert2015) for systems where heat fluxes are negligible and are well magnetised so the particle gyroradius is much smaller than the relevant spatial gradients present in the plasma. This provides an interesting framework that can be relevant for general relativistic MHD simulations used in astrophysical studies of, e.g. low-luminosity accretion disks around SMBH (EHT Collaboration 2019), pressure anisotropy effects on cosmic-ray hydrodynamics (Zweibel Reference Zweibel2020) and for extending CGL-MHD fluid models like those used for studying magneto-immutability effects (Majeski et al. Reference Majeski, Kunz and Squire2024) in the ICM of galaxy clusters, in order to include relativistic species that are part of the ICM (e.g. cosmic rays). Additionally, the linear stability of our time-dependent anisotropic Maxwell–Jüttner distribution can then be studied analytically and numerically (e.g. using linear solvers like the Arbitrary Linear Plasma Solver (Verscharen et al. Reference Verscharen, Klein, Chandran, Stevens, Salem and Bale2018)), along with the onset and evolution of pressure anisotropy-driven microinstabilities relevant in the relativistic regime. The results can then be compared directly with driven PIC simulations. A direct implementation of this anisotropic Maxwell–Jüttner distribution in PIC simulations is also possible, and will be addressed in a future work.
Finally, our general solution to the drift kinetic equation also admits other types of interesting initial distributions, including power-law and kappa distributions, relevant for cosmic-ray and solar wind studies, respectively. Double-adiabatic equations and linear stability analysis can also be performed in each case.
Acknowledgements
The authors thank the anonymous referees for the insightful comments and suggestions that improved the paper. We thank M. Riquelme, A. Sandoval and M. Astete for useful comments and discussions. During the final stages of preparation of this manuscript, we became aware of a similar and independent study by A. Wierzchucka, P. J. Bilbao, A. G. R. Thomas, D. A. Uzdensky and A. A. Schekochihin (accepted for publication in JPP), carried out independently and simultaneously to ours (Wierzchucka et al. Reference Wierzchucka, Bilbao, Thomas, Uzdensky and Schekochihin2026). The authors thank Vedant Dhruv for pointing out a typo in the exponent in one of the equations, that was corrected in this version of the paper.
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
F.L acknowledges support by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project 24-612 (‘Excitation and dissipation of kinetic-scale fluctuations in space plasmas’). A.T. was supported by the Department of Energy Fusion Energy Sciences Postdoctoral Research Program, administered by the Oak Ridge Institute for Science and Education and Oak Ridge Associated Universities under DOE contract DE-SC0014664. E.G.Z. acknowledges support from the National Science Foundation through Grant PHY-2409224 and the Simons Foundation through the Simons Collaboration on Extreme Electrodynamics of Compact Sources as well as the hospitality of the Leonard Parker Center for Gravitation, Cosmology, and Astrophysics at UW–Milwaukee, where some of this work was completed.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Mass ratio and plasma beta dependence in our PIC simulations
A.1. The
$m_i/m_e$
and
$\beta _i^{\text{init}}$
dependency in shearing simulations
In our shearing simulations, we tested various cases for both the ion to electron mass ratio
$m_i/m_e$
and initial plasma beta
$\beta _i^{\text{init}}$
. We saw no dependency on any of these two parameters, analogous to the non-relativistic CGL case. These results are collected in figure 7. We are showing ion quantities only, but the evolution of electron pressures is analogous.
For the mass ratio, we tried three different values:
$m_i/m_e = 8, 64\, {\rm and}\,1836$
(panels
$a$
and
$b$
in figure 7). For both ion perpendicular and parallel pressures, we can see that all three mass ratio cases are well described by their respective ultrarelativistic evolution (cf. equations (2.18), (2.21)).
Similarly, for the initial beta, we tried three different values spanning two orders of magnitude:
$\beta _i^{\text{init}}= 0.05, 0.5$
and
$ 5$
(panels
$c$
and
$d$
in figure 7). For both ion perpendicular and parallel pressures, their evolution is well explained by their respective ultrarelativistic evolution in all three cases.
(
$a$
) The evolution of the ion perpendicular pressure for runs Shb0.5m8d30wci3200 (pale blue circles), Shb0.5m64d30wci3200 (orange triangles) and Shb0.5m1836d30wci3200 (green stars). (
$b$
) The evolution of the ion parallel pressure for runs Shb0.5m8d30wci3200 (pale blue circles), Shb0.5m64d30wci3200 (orange triangles) and Shb0.5m1836d30wci3200 (green stars). (
$c$
) The evolution of the ion perpendicular pressure for runs Shb0.05m8d30wci3200 (blue circles), Shb0.5m8d30wci3200 (red triangles) and Shb5m8d30wci3200 (green stars). (
$d$
) The evolution of the ion parallel pressure for runs Shb0.05m8d30wci3200 (blue circles), Shb0.5m8d30wci3200 (red triangles) and Shb5m8d30wci3200 (green stars). In panels
$(a)$
and
$(c)$
, the non-relativistic CGL evolution for
$P_{\perp }$
is shown in dashed grey line, and the ultrarelativistic evolution for
$P_{\perp }$
(cf. equation (2.18)) is shown in solid black line. In panels
$(b)$
and
$(d)$
, the non-relativistic CGL evolution for
$P_{\parallel }$
is shown in dashed grey line and the ultrarelativistic evolution for
$P_{\parallel }$
(cf. equation (2.21)) is shown in solid black line.
Figure 7. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the solid lines fully coincide with the open circles, open triangles, and open stars, and the gray dashed line agree only at the beginning, and diverges at later times. Read caption for more details.
A.2. The
$m_i/m_e$
and
$\beta _i^{\text{init}}$
dependency in compressing simulations
For our compressing simulations, we also tested the dependency of our results with mass ratio and initial plasma beta. We tried the same values as in the shearing case, namely,
$m_i/m_e = 8, 64$
and
$1836$
and
$\beta _i^{\text{init}} = 0.05, 0.5$
and
$5$
. All simulations have initial
$k_BT/m_ic^2=30$
, and we are only showing ion results, as electron show a analogous behaviour. The results are shown in figure 8.
(
$a$
) The evolution of the ion perpendicular temperature for runs Compb0.5m8d30wci3200 (pale blue circles), Compb0.5m64d30wci3200 (orange triangles) and Compb0.5m1836d30wci3200 (green stars). (
$b$
) The evolution of the ion parallel temperature for runs Compb0.5m8d30wci3200 (pale blue circles), Compb0.5m64d30wci3200 (orange triangles) and Compb0.5m1836d30wci3200 (green stars). (
$c$
) The evolution of the ion perpendicular temperature for runs Compb0.05m8d30wci3200 (blue circles), Compb0.5m8d30wci3200 (red triangles) and Compb5m8d30wci3200 (green stars). (
$d$
) The evolution of the ion parallel temperature for runs Compb0.05m8d30wci3200 (blue circles), Compb0.5m8d30wci3200 (red triangles) and Compb5m8d30wci3200 (green stars). In panel
$(a)$
and
$(c)$
, the non-relativistic CGL evolution for
$T_{\perp }$
is shown in dashed grey line, and the ultrarelativistic evolution for
$T_{\perp }$
(cf. (4.4)) is shown in solid black line. In panels
$(b)$
and
$(d)$
, the non-relativistic CGL evolution for
$T_{\parallel }$
is shown in dashed grey line, and the ultrarelativistic evolution for
$T_{\parallel }$
(cf. (4.6)) is shown in solid black line.
Figure 8. Long description
A mosaic of four graphs describing the evolution of physical quantities described in the caption for a specific simulation run. In each of the four panels, the solid lines fully coincide with the open circles, open triangles, and open stars, and the gray dashed line does not agree with the rest of the curves at all. Read caption for more details.
We can see that for the three mass ratio considered, the evolutions of the perpendicular and parallel pressures (panels
$a$
and
$b$
in figure 8) are well described by the ultrarelativistic (4.4) and (4.6), respectively. We observe a similar pattern for the simulations with different initial beta (panels
$c$
and
$d$
in figure 8). In this case, both
$T_{\perp ,i}$
and
$T_{\parallel ,i}$
evolutions are well explained by (4.4) and (4.6) for all the values of beta considered. Therefore, similar to the non-relativistic case, the relativistic extension of the double-adiabatic equation shows a robust independence of both mass ratio and initial plasma beta. This independence of plasma beta allows the application of these extensions to the double-adiabatic equations to a wide variety of astrophysical scenarios in which one or more species are relativistically hot.
βiinit=8πnikBTiinit/B2
mi/me
kBTiinit/mic2
Ωi/s
Nppc=300
d/Δx≡c/(ωp,e2+ωp,i2)1/2/Δx
L/RL,i
βiinit=8πnikBTiinit/B2
mi/me
kBTiinit/mic2
Ωi/q
Nppc=300
d/Δx≡c/(ωp,e2+ωp,i2)1/2/Δx
L/RL,i
a
kBTiinit/mic2=0.2
P⊥
b
P∥
c
kBTeinit/mec2=1.6
mi/me=8
P⊥
d
P∥
a
c
P⊥
b
d
P∥
a
kBTiinit/mic2=30
P⊥
P⊥
b
P∥
P∥
c
kBTeinit/mec2=55080
P⊥
P⊥
d
P∥
P∥
(a)
(c)
P⊥
(b)
(d)
P∥
a
kBTiinit/mic2=0.2
T⊥
b
T∥
c
kBTeinit/mec2=1.6
T⊥
(d)
T∥
(a)
(c)
P⊥
(b)
(d)
P∥
a
kBTiinit/mic2=30
T⊥
b
T∥
(a)
(c)
P⊥
(b)
(d)
P∥
a
t⋅s=0
t⋅s=1
faMJ(p⊥,t)
t⋅s=0
t⋅s=1
θ=0.2
b
t⋅s=0
t⋅s=1
faMJ(p∥,t)
t⋅s=0
t⋅s=1
θ=0.2
(c)
(d)
(a)
(b)
θ=30
a
t⋅s=0
t⋅s=1
faMJ,c(p⊥,t)
t⋅s=0
t⋅s=1
θ=0.2
b
t⋅s=0
t⋅s=1
faMJ,c(p⊥,t)
t⋅s=0
t⋅s=1
θ=30
a
b
c
d
(a)
(c)
P⊥
P⊥
(b)
(d)
P∥
P∥
a
b
c
d
(a)
(c)
T⊥
T⊥
(b)
(d)
T∥
T∥
Open access
