1. Introduction
In magnetised, collisionless plasmas, constituent particles gyrate rapidly around magnetic-field lines while experiencing infrequent Coulomb collisions. Consequently, the particle gyration frequency sets the shortest time scale in the system, whereas the Coulomb collision time is the largest. Studies of collisionless and weakly collisional magnetised plasmas in the non-relativistic regime have provided significant insight into the dynamics governing space and astrophysical environments such as the solar wind (Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009), the Earth’s magnetosphere (Øieroset et al. Reference Øieroset, Phan, Fujimoto, Lin and Lepping2001; Burch et al. Reference Burch, Moore, Torbert and Giles2016), and the intracluster medium (Schekochihin & Cowley Reference Schekochihin and Cowley2006; Kunz, Jones & Zhuravleva Reference Kunz, Jones and Zhuravleva2022). Collisionless plasmas are also abundant in high-energy relativistic astrophysical systems, like pulsar-wind nebulae (Kargaltsev et al. Reference Kargaltsev, Cerutti, Lyubarsky and Striani2015), pulsar magnetospheres (Cerutti & Beloborodov Reference Cerutti and Beloborodov2017; Philippov & Kramer Reference Philippov and Kramer2022), accretion flows onto black holes (Rees et al. Reference Rees, Begelman, Blandford and Phinney1982; Narayan & Yi Reference Narayan and Yi1995; Quataert Reference Quataert2003; Yuan & Narayan Reference Yuan and Narayan2014), and jets from active galactic nuclei (e.g. blazars) (Begelman, Blandford & Rees Reference Begelman, Blandford and Rees1984; Blandford, Meier & Readhead Reference Blandford, Meier and Readhead2019). Many of these sources exhibit pronounced non-thermal emission spectra, providing direct observational evidence that the associated plasmas are far from Maxwellian equilibrium. Such departure from thermal equilibrium means that studying collisionless magnetised plasmas is challenging and requires a kinetic approach.
Rapid gyromotion renders the particle distribution function azimuthally symmetric in velocity space about the magnetic field, a property known as gyrotropy. The pressure tensor in such a plasma is diagonal. When the magnetic field varies on time scales much longer than the gyration period, the first adiabatic invariant,
$\mu$
, which is the angular momentum of a particle on a Larmor orbit, is conserved (Gardner Reference Gardner1959), leading to a non-identical evolution of the perpendicular and parallel pressures,
$P_\perp$
and
$P_\parallel$
. In the non-relativistic limit, this evolution is described by the Chew–Goldberger–Low (CGL) equations (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956), which express the Lagrangian evolution of
$P_\perp$
and
$P_\parallel$
in terms of the particle number density,
$n$
, and magnetic-field strength,
$B$
, as well as heat fluxes. When the latter can be neglected, the CGL equations are known as the double-adiabatic approximation and provide a practical fluid closure of the system: in a moving fluid element,
$P_\perp$
and
$P_\parallel$
are linked to the density and magnetic-field strength by instantaneous adiabatic equations of state (EoSs),
$P_\perp \propto nB$
and
$P_\parallel \propto n^3/B^2$
. These relations are the anisotropic generalisation of the familiar adiabatic EoS
$P \propto n^{\varGamma }$
of a highly collisional, isotropic plasma, where the adiabatic index
$\varGamma$
ranges from
$5/3$
in the non-relativistic regime to
$4/3$
in the ultra-relativistic limit (Synge Reference Synge1957).
When heat fluxes can be neglected, the CGL equations are sometimes referred to as ‘double adiabatic’ because they encompass the conservation of two adiabatic invariants, the first being the aforementioned
$\mu$
. The second adiabatic invariant,
$J$
, is associated with the trapping and bouncing of particles in magnetic mirrors and remains conserved under the condition that this bouncing is much faster than the fluid motions. Such a requirement is much more restrictive than that of rapid gyration, and so the conservation of
$J$
is less robust than the conservation of
$\mu$
. The conservation of these two quantities can be shown through the use of action-angle variables (Gardner Reference Gardner1959; Kruskal Reference Kruskal1962) or directly by expanding the particle equations of motion in small gyration periods, an approximation known as drift kinetics (Northrop & Teller Reference Northrop and Teller1960; Vandervoort Reference Vandervoort1960; Northrop Reference Northrop1963).
The double-adiabatic equations give insight into the various processes governing the dynamics of collisionless plasmas. As they predict different evolution of
$P_\perp$
and
$P_\parallel$
, local changes in the density or magnetic field, e.g. due to shearing or compressing flows, can give rise to local pressure anisotropies (Hellinger & Trávníček Reference Hellinger and Trávníček2008; Kunz, Schekochihin & Stone Reference Kunz, Schekochihin and Stone2014; Sironi & Narayan Reference Sironi and Narayan2015). In high-
$\beta$
plasmas, where
$\beta$
is the ratio of the thermal to magnetic pressure, pressure anisotropy can lead to the excitation of various kinetic instabilities (see Bott, Cowley & Schekochihin Reference Bott, Cowley and Schekochihin2024 and references therein). When the anisotropy,
$\varDelta \equiv P_\perp /P_\parallel - 1$
, is positive, the mirror (Hasegawa Reference Hasegawa1969; Southwood & Kivelson Reference Southwood and Kivelson1993) and ion-cyclotron instabilities emerge if
$\varDelta \gtrsim 1/\beta$
. In the opposite case, if
$\varDelta \lesssim -2/\beta$
, the microscale dynamics of the plasma becomes dominated by the firehose instability (Chandrasekhar, Kaufman & Watson Reference Chandrasekhar, Kaufman and Watson1958). After a short period of exponential growth, both the mirror and firehose instabilities saturate and, albeit through different mechanisms, pin the pressure anisotropy to their respective thresholds (Melville, Schekochihin & Kunz Reference Melville, Schekochihin and Kunz2016; Bott et al. Reference Bott, Kunz, Quataert, Squire and Arzamasskiy2025). As a result, the pressure becomes nearly isotropic and evolves with an adiabatic index
$\varGamma = 5/3$
in the high-
$\beta$
non-relativistic regime.
As the distribution function is gyrotropic to lowest order in the gyration period, the Vlasov equation can be averaged over the gyromotion to obtain the so-called drift-kinetic equation. This constitutes a starting point for studying the kinetic instabilities excited in collisionless plasmas, although finite-gyroradius effects must be included to recover the fastest-growing modes (Yoon, Wu & de Assis Reference Yoon, Wu and de Assis1993; Hellinger Reference Hellinger2007; Bott et al. Reference Bott, Cowley and Schekochihin2024).
The CGL equations are derived by taking perpendicular and parallel pressure moments of the drift-kinetic equation (see e.g. Schekochihin et al. Reference Schekochihin, Cowley, Rincon and Rosin2010 for a derivation). This procedure can be extended to relativistic collisionless plasmas (Gedalin & Oiberman Reference Gedalin and Oiberman1995; Hegade & Stone Reference Hegade and Stone2026; Wierzchucka et al. Reference Wierzchucka, Bilbao, Uzdensky and Schekochihin2026b
). Unlike in the non-relativistic regime, the relativistic CGL equations do not provide an evolution equation for the two pressures. Instead, they evolve related quantities that we refer to below as the modified pressures. However, as it is the pressures, and not the modified pressures, that are needed in the fluid momentum equation, this means that, relativistically, the generalised CGL equations cannot provide a sufficient fluid closure of the system. Describing the evolution of
$P_\perp$
and
$P_\parallel$
in a relativistic, collisionless, magnetised plasma is somewhat non-trivial. The main reason for this is that, in the relativistic case, the momentum-space integration in the definition of
$P_\perp$
and
$P_\parallel$
involves, in addition to the corresponding momentum components, the particle Lorentz factor
$\gamma$
, which is a nonlinear function of the perpendicular and parallel momenta. As a result,
$P_\perp$
(
$P_\parallel$
) depends not only on the perpendicular (parallel) component of the particle momentum but also on the total momentum (through the Lorentz factor). In other words, each pressure component depends on the full momentum distribution rather than solely on its corresponding momentum projection, as in the non-relativistic case. Therefore, their evolution in the double-adiabatic limit, i.e. in the absence of heat fluxes, turns out to depend not only on the density and magnetic-field strength but also on the pressure anisotropy.
In this paper, we derive the double-adiabatic equations satisfied by
$P_\perp$
and
$P_\parallel$
in relativistic collisionless plasmas, with a non-relativistic bulk flow. We do this by abandoning the standard moment-based derivation of the evolution equations and instead formulating a theory of adiabatic invariance focused on the system’s symmetries. Through the use of phase-space volume conservation properties, we show that imposing certain symmetries leads to a self-similar evolution of the distribution function. Calculating the pressure moments of the evolved distribution then yields the adiabatic EoSs, whose form, in general, depends on the distribution at some initial (reference) time. For example, one can easily show that a distribution function isotropic in peculiar momentum satisfies
$P \propto n^{5/3}$
and
$P \propto n^{4/3}$
in the non-relativistic and ultra-relativistic regimes, respectively. If the distribution function is gyrotropic and has parity symmetry along the magnetic-field direction in the peculiar-momentum coordinate, we find by the same procedure that the double-adiabatic equations,
$P_\perp \propto nB$
and
$P_\parallel \propto n^3/B^2$
, hold in the non-relativistic limit.
Using this symmetry-based approach, we extend the double-adiabatic equations to the relativistic regime. We study these equations for the case of an ultra-relativistically hot plasma with an isotropic reference distribution. We find that the exact form of the double-adiabatic equations depends on the pressure anisotropy and derive the asymptotic theory for three main limiting cases. When
$P_\perp \gg P_\parallel$
(
$\varDelta \gg 1)$
, we find
$P_\perp \propto n B^{1/2}$
and
$P_\parallel \propto n^3/B^{5/2}$
, whereas when
$P_\parallel \gg P_\perp$
(
$\varDelta \approx -1$
), we find
$P_\perp \propto B^2 \ln (n^{\prime 2}/B^{\prime 3})$
and
$P_\parallel \propto {n^2/B}$
. Here,
$B'$
and
$n'$
are, respectively, the relative changes in the magnetic-field strength and density. These evolution equations agree with the initial discussion of the relativistic CGL equations in Gedalin & Oiberman (Reference Gedalin and Oiberman1995), up to the logarithmic correction in the second
$P_\perp$
expression. In contrast, close to isotropy, i.e. when
$|\varDelta | \ll 1$
, we find
$P_\perp \propto (nB)^{4/5}$
and
$P_\parallel \propto (n^{3}/B^{2})^{4/5}$
. Our results differ drastically from the non-relativistic theory (Chew et al. Reference Chew, Goldberger and Low1956), in which the same double-adiabatic relations hold for all anisotropies. We confirm our predictions using two-dimensional particle-in-cell (PIC) simulations incorporating a large-scale compressive flow.
In situations where heat fluxes are negligible, these double-adiabatic equations can be the EoSs for relativistic collisionless plasmas. Such EoSs are crucial for modelling large-scale high-energy astrophysical systems by means of a fluid theory, e.g. to determine regions where kinetic processes like the firehose and mirror instabilities are excited in high-
$\beta$
plasmas. Their applicability to the relativistic regime also means that our EoSs can be coupled to other high-energy processes, such as synchrotron cooling, providing a pathway to modelling phenomena like the recently discovered synchrotron firehose instability (Zhdankin, Kunz & Uzdensky Reference Zhdankin, Kunz and Uzdensky2023). Double-adiabatic EoSs may additionally have an important effect on the large-scale magnetohydrodynamic-equilibrium structure of astrophysical magnetospheres and flows (Beskin & Kuznetsova Reference Beskin and Kuznetsova2000).
2. Adiabatic invariance and symmetries
We consider a particle distribution function
$f(\boldsymbol{r}, \boldsymbol{p}, t)$
, where
$\boldsymbol{r}$
and
$\boldsymbol{p}$
denote the position and momentum variables, and
$t$
denotes time. The distribution function is the density of a phase-space fluid filling the six-dimensional
$(\boldsymbol{r}, \boldsymbol{p})$
space, moving with velocity
The velocity coordinate
$\boldsymbol{v}$
is related to the momentum through
$\boldsymbol{p} = \gamma m \boldsymbol{v}$
, where
$\gamma = \sqrt {1+(p/mc)^2}$
is the Lorentz factor, and
$m$
is the particle mass. The force
$\boldsymbol{F}$
depends on the particular physical system that the distribution function describes. In this paper, we focus on Hamiltonian systems, in which case the phase-space fluid has the additional property of being incompressible.
To describe how
$f$
evolves in time, we adopt a Lagrangian formalism (Newcomb Reference Newcomb1961). We describe the position of an infinitesimal fluid element in phase space using not its current coordinates
$(\boldsymbol{r}, \boldsymbol{p})$
, but instead its position at
$t=0$
,
$(\boldsymbol{r}_0, \boldsymbol{p}_0)$
. In this approach, the conservation of probability is
and the incompressibility of the phase-space fluid, also referred to as the conservation of phase-space volume, is written as
One can combine (2.2) and (2.3) to find Liouville’s theorem
where
$f_0$
is the distribution function at
$t=0$
.
Two important moments of
$f$
are the density and the Eckart (Reference Eckart1940) bulk velocity, defined, respectively, by
In this work, for simplicity, we restrict our analysis to the non-relativistic limit as far as the bulk velocity of the plasma is concerned,
$|\boldsymbol{u}| \ll c$
. Nevertheless, our findings regarding the adiabatic and double-adiabatic EoSs remain valid in the fully relativistic regime, as discussed in § 5. We label quantities in the comoving frame with a tilde, referring to them as ‘peculiar’. To first order in
$|\boldsymbol{u}|/c$
, the Lorentz transformation of the momentum to the peculiar momentum is
and the corresponding transformation of the Lorentz factor is
The Jacobian of this change of variables is, again to first order in
$|\boldsymbol{u}|/c$
,
where
$\tilde {\boldsymbol{v}} = \tilde {\boldsymbol{p}}/\tilde {\gamma }m$
. In the fully non-relativistic regime,
$|\tilde {\boldsymbol{p}}| \ll m c$
, (2.6) reduces to
$\boldsymbol{p} = m\boldsymbol{u} + \tilde {\boldsymbol{p}}$
, as expected. Note that, while by definition of the bulk velocity (2.5),
$\int {\rm d}\tilde {\boldsymbol{p}} \ \tilde {\boldsymbol{v}}f = 0$
, in general
$\int {\rm d}\tilde {\boldsymbol{p}} \ \tilde {\boldsymbol{p}}f \neq 0$
, reflecting the differences between the Landau (Landau & Lifshitz Reference Landau and Lifshitz1987) and Eckart bulk velocities.
2.1. Isotropy
For illustration, we first consider a distribution function that evolves in such a way as to remain isotropic in the frame moving with
$\boldsymbol{u}$
. This could arise, e.g. due to efficient pitch-angle scattering dominating the system’s dynamics, and implies the neglect of heat fluxes. Introducing spherical polar coordinates
$(\tilde {p}, \theta , \phi )$
, with the polar axis orientated in an arbitrary direction, we have
and the corresponding Jacobian is
The isotropy assumption is then the statement that
We assume that in this ‘reduced phase space’, phase-space volume is conserved, i.e.
where we have made use of (2.10). Combining (2.12) with Liouville’s theorem (2.4) then implies the conservation of probability
Considering a rapidly isotropised distribution function is equivalent to assuming that the phase space can be described by its average over some fast time scale
$t_{\text{fast}}$
associated with the isotropisation process. In the laboratory frame, such an average is defined by
\begin{align} \langle \boldsymbol{\cdot }\rangle _{t_{\textrm {fast}}} \equiv \frac {\int _t^{t+t_{\textrm {fast}}}(\boldsymbol{\cdot })\,{\rm d}t} {\int _t^{t+t_{\textrm {fast}}}{\rm d}t}. \end{align}
However, isotropisation is a microscopic process that occurs in the particle’s proper time
$\tau$
, which is related to the laboratory time via
${\rm d}t=\gamma \,{\rm d}\tau$
. We assume that, during this process, the particle’s motion becomes uniformly randomised as far as its direction is concerned, allowing the proper-time average to be replaced by an average over the solid angle of the peculiar momentum. Therefore,
\begin{align} \langle \boldsymbol{\cdot }\rangle _{t_{\textrm {fast}}} = \frac {\int _\tau ^{\tau +\tau _{\textrm {fast}}}(\boldsymbol{\cdot })\gamma \,{\rm d}\tau } {\int _\tau ^{\tau +\tau _{\textrm {fast}}}\gamma \,{\rm d}\tau } = \frac {\iint (\boldsymbol{\cdot })\gamma \,\sin \theta \, {\rm d} \theta {\rm d}\phi } {\iint \gamma \,\sin \theta \, {\rm d} \theta {\rm d}\phi } = \frac {\langle \gamma (\boldsymbol{\cdot })\rangle _{\theta ,\phi }} {\langle \gamma \rangle _{\theta ,\phi }}. \end{align}
Thus, in the relativistic limit, replacing the averaging over fast time scales with angle averaging leads to the appearance of a Lorentz factor inside the angle integrals. As expected, in the non-relativistic limit, where
$\gamma \approx 1$
, these averages are the same.
Averaging the phase-space flow (2.1) via (2.15) gives
$\langle \dot {\boldsymbol{r}}\rangle _{t_{\text{fast}}} = \boldsymbol{u}$
, where we have used (2.6) and (2.7). Consequently, in the reduced phase space, the mapping between
$\boldsymbol{r}$
and
$\boldsymbol{r}_0$
is independent of the peculiar coordinates. The phase-space path
$\boldsymbol{r} = \boldsymbol{r}(t, \boldsymbol{r}_0)$
then corresponds to that of a fluid element advected with the bulk velocity
$\boldsymbol{u}$
in three-dimensional position space. The Jacobian (2.12) therefore factorises into spatial and momentum parts
We can calculate the first factor in (2.16) by returning to (2.13), integrating over the momenta, and recalling the definition of density (2.5) to find
which implies
Combining the above expression with (2.16) and inserting the result into (2.12) yields a differential equation whose solution is
We can now substitute (2.19) into Liouville’s theorem (2.4) and find the evolution of the distribution function
Note that the explicit form of
$\boldsymbol{r}(t, \boldsymbol{r}_0)$
is not required, as the trajectory corresponds to that of a fluid element advected with the bulk velocity
$\boldsymbol{u}$
. Hence, any evolution equations obtained by taking moments of (2.20) are Lagrangian relations with respect to
$\boldsymbol{u}$
.
2.1.1. Equation of state for isotropic plasma
Now that we have found how the distribution function evolves, we can work out the corresponding EoS. In the relativistic regime, pressure is given by
Inserting (2.20) into (2.21), integrating over
$\theta$
and
$\phi$
, and changing the momentum integration variable to
$\tilde {p}_0$
via (2.19) gives the adiabatic EoS for an isotropic plasma
\begin{align} P = \frac {4 \pi n^{\prime {5/3}}}{3m} \int {\rm d}\tilde {p}_0 \, \frac {\tilde {p}^4_0 f_0(\tilde {p}_0)} {\sqrt {1 + n^{\prime {2/3}} \left ( {\tilde {p}_0}/{mc} \right )^2 }}. \end{align}
In the non-relativistic regime,
$\tilde {p} \ll mc$
, (2.22) reduces to the familiar scaling
$P \propto n^{5/3}$
and in the ultra-relativistic limit,
$\tilde {p} \gg mc$
, we instead obtain
$P \propto n^{4/3}$
, as expected. In general, however,
$P$
depends not solely on density but also on the form of the initial distribution function
$f_0(\tilde {p})$
. Substituting the Maxwell–Jüttner distribution into (2.22) recovers the Synge gas EoS (Synge Reference Synge1957).
2.2. Gyrotropy and parity
We now consider a more general case, where the isotropy assumption is relaxed but the distribution function still satisfies two symmetries in the peculiar-momentum coordinate. First, it is gyrotropic, i.e. it has cylindrical symmetry with respect to some vector field, which we call
$\boldsymbol{B}$
. Second, the distribution function has parity symmetry along this vector. We assume that the vector
$\boldsymbol{B}$
is divergence free and frozen into the bulk flow; it therefore satisfies the induction equation
For example, in a drift-kinetic plasma,
$\boldsymbol{B}$
corresponds to the magnetic field (Kulsrud Reference Kulsrud1964). A distribution function with such symmetries can arise in systems where fast processes, such as gyromotion or the bouncing of trapped particles in magnetic mirrors of collisionless plasmas, dominate the dynamics, and implies the neglect of heat fluxes.Footnote
1
We adopt a cylindrical coordinate system in the peculiar-momentum space
$(\tilde {p}_\perp , \phi , \tilde {p}_\parallel , \sigma )$
, with the polar axis parallel to
$\boldsymbol{b} \equiv \boldsymbol{B}/B$
where
$\tilde {p}_\parallel \geqslant 0$
and
$\sigma$
takes the values
$+1$
or
$-1$
, labelling the direction of the parallel peculiar momentum. The Jacobian of this change of variables is
The gyrotropy and parity assumptions are then equivalent to imposing that
$f$
is independent of
$\phi$
and
$\sigma$
We again assume that phase-space volume in the reduced phase space is conserved, viz.
where we have made use of (2.25). Combining Liouville’s theorem (2.4) with (2.27) yields the conservation of probability (2.2) in our reduced phase space
In analogy to § 2.1, imposing (2.26) is equivalent to assuming that the phase space is well described by its average over some fast time scales associated with the two symmetries (e.g. Larmor gyration and bouncing inside magnetic mirrors). In the non-relativistic limit, averaging over these time scales in the laboratory frame is equivalent to an average over
$\sigma$
and
$\phi$
. In the relativistic case, as in (2.15), the average acquires a weight of
$\gamma$
\begin{align} \langle \boldsymbol{\cdot }\rangle _{t_{\text{fast}}} \equiv \frac {\int _{t}^{t + t_{\text{fast}}} (\boldsymbol{\cdot }) {\rm d}t}{\int _{t}^{t + t_{\text{fast}}}{\rm d}t} = \frac {\sum _\sigma \int (\boldsymbol{\cdot })\gamma \,{\rm d}\phi } {\sum _\sigma \int \gamma \,{\rm d}\phi } = \frac {\langle \gamma (\boldsymbol{\cdot })\rangle _{\phi , \sigma }} {\langle \gamma \rangle _{\phi ,\sigma }}. \end{align}
Averaging the phase-space flow (2.1) over the fast time scales using (2.6), (2.7), and (2.29) gives
$\langle \dot {\boldsymbol{r}}\rangle _{t_{\text{fast}}} = \boldsymbol{u}$
, as in § 2.1, implying that, in the reduced phase space, the trajectory
$\boldsymbol{r} = \boldsymbol{r}(t, \boldsymbol{r}_0)$
again corresponds to that of a fluid element advected with the bulk flow
$\boldsymbol{u}$
. Hence, the Jacobian in (2.27) factorises into spatial and momentum contributions
where the former can be evaluated using (2.18). Then, from (2.30) and (2.27), we have
To determine the evolution of
$f$
, we require expressions for
$\tilde {p}_\perp$
and
$\tilde {p}_\parallel$
in terms of their initial values. However, we have obtained only a single constraint (2.31). We therefore assume that, in addition to phase-space volume, the phase-space area perpendicular to
$\boldsymbol{B}$
is conserved. A natural coordinate system for resolving directions parallel and perpendicular to
$\boldsymbol{B}$
is provided by the Clebsch variables
$(\alpha , \beta , \ell )$
, defined by
Here,
$\alpha$
and
$\beta$
label a magnetic-field line, while
$\ell$
measures distance along it. The Jacobian of this transformation is
and the conservation of perpendicular phase-space area can be written as
or, equivalently,
As we discuss in § 2.2.3, this is just the conservation of
$\mu$
in collisionless, magnetised plasmas, or rather the conservation of
$\mu$
is how the conservation of perpendicular phase-space area is physically achieved by a plasma.
From the induction equation (2.23), we have
$\dot {\alpha } = \dot {\beta } = 0$
, so
$\alpha = \alpha _0$
and
$\beta = \beta _0$
, i.e.
$\alpha$
and
$\beta$
remain constant as the field is advected by the flow
$\boldsymbol{u}$
. Hence, (2.35) yields
Integration of (2.36) yields
and substituting this into (2.31) gives
Finally, we can combine (2.37) and (2.38) with Liouville’s theorem (2.4) to find the evolved distribution function
Again, the time dependence is encoded via
$n'(t, \boldsymbol{r}_0)$
,
$B'(t, \boldsymbol{r}_0)$
, and
$\boldsymbol{r}(t, \boldsymbol{r}_0)$
.
2.2.1. Equations of state for gyrotropic plasma
Similarly to § 2.1.1, we can now use the evolved distribution (2.39) to calculate the corresponding EoSs. For a relativistic plasma, the parallel and perpendicular pressures are defined by
Inserting (2.39) into these equations, integrating over
$\sigma$
and
$\phi$
, and finally changing the integration variables to
$\tilde {p}_{\perp 0}$
and
$\tilde {p}_{\parallel 0}$
via (2.37) and (2.38) yields the double-adiabatic pressure-evolution equations
\begin{align} P_\perp &= \frac {{2\pi } n' B'}{m}\int {\rm d}\tilde {p}_{\perp 0} \ {\rm d}\tilde {p}_{\parallel 0} \ \nonumber \\[5pt] & \frac {\tilde {p}_{\perp 0}^3 f_0(\tilde {p}_{\perp 0}, \tilde {p}_{\parallel 0}) \ } {\sqrt {1 + B' \left ({\tilde {p}_{\perp 0}}/{mc} \right )^2 + \left ( {n'}/{B'}\right )^2({\tilde {p}_{\parallel 0}}/{mc} )^2}},\\[-10pt]\nonumber \end{align}
\begin{align} P_\parallel &= \frac {{4\pi }n^{\prime 3}}{m B^{\prime 2}}\int {\rm d}\tilde {p}_{\perp 0} \ {\rm d}\tilde {p}_{\parallel 0} \ \nonumber \\[5pt] & \frac {\tilde {p}_{\perp 0} \tilde {p}_{\parallel 0}^2 f_0(\tilde {p}_{\perp 0}, \tilde {p}_{\parallel 0}) \ } {\sqrt {1 + B' \left ({\tilde {p}_{\perp 0}}/{mc} \right )^2 + \left ( {n'}/{B'}\right )^2({\tilde {p}_{\parallel 0}}/{mc} )^2}}. \end{align}
These integrals do not, in general, exhibit a single set of simple scalings with
$n'$
and
$B'$
for arbitrary
$f_0$
. However, in the non-relativistic limit,
$\tilde {p}_\perp , \tilde {p}_\parallel \ll mc$
, (2.41) and (2.42) reduce to the CGL equations, i.e.
$P_\perp \propto n B$
and
$P_\parallel \propto n^3/B^2$
(Chew et al. Reference Chew, Goldberger and Low1956), for all
$f_0$
.
2.2.2. Modified pressures
The complicated dependence of (2.41) and (2.42) on
$n'$
and
$B'$
arises due to the presence of the Lorentz factor in the integrand of (2.40). To avoid this issue, one can define the modified pressures
Proceeding as before and inserting (2.39) into the definitions above, we find
The integrals in both (2.44) and (2.45) can be replaced by the initial values of the modified pressures,
$\hat {P}_{\perp 0}$
and
$\hat {P}_{\parallel 0}$
. We thus obtain the simple relationships
where
$\hat {P}_\perp ' \equiv \hat {P}_{\perp }/\hat {P}_{\perp 0}$
and
$\hat {P}_\parallel ' \equiv \hat {P}_{\parallel }/\hat {P}_{\parallel 0}$
. These expressions are the relativistic generalisation of the CGL equations. They can alternatively be derived by taking moments of the relativistic Vlasov equation and neglecting heat fluxes, as was first done in Gedalin & Oiberman (Reference Gedalin and Oiberman1995). However, in the fluid momentum equation, it is
$P_\perp$
and
$P_\parallel$
that appear, not
$\hat {P}_\perp$
and
$\hat {P}_\parallel$
. Hence, in general, (2.46) does not provide a suitable fluid closure for collisionless relativistic plasmas.
2.2.3. Application to collisionless, magnetised plasmas
The theory presented above derives double-adiabatic evolution equations for a fluid system with a divergence-free vector field
$\boldsymbol{B}$
frozen into the bulk flow. The derivation relies on a set of assumptions regarding the symmetries of the distribution function in peculiar-momentum space, as well as the conservation properties of the evolving phase-space density. First, we considered a gyrotropic distribution function and assumed conservation of phase-space area perpendicular to
$\boldsymbol{B}$
. Second, we imposed parity symmetry along the direction of
$\boldsymbol{B}$
in peculiar-momentum space, together with conservation of the full phase-space volume.
An example of a system in which such symmetries arise as a result of fast processes is a collisionless magnetised plasma. The fast particle gyration provides the largest frequency in such a system, implying the gyrotropy of the distribution function and the conservation of the first adiabatic invariant (Gardner Reference Gardner1959)
along gyro-averaged particle trajectories. This conservation implies (2.37) and provides the physical mechanism whereby the conservation of perpendicular phase-space area (2.34) holds. We are currently not aware of a system where such conservation is achieved by different means.
The conservation of
$\mu$
can cause particles with sufficiently large pitch angles to become trapped in local minima of the magnetic-field strength (magnetic mirrors) and undergo bounce motion. If such bouncing occurs on time scales much shorter than those of the fluid motion, a second adiabatic invariant,
is conserved, where the integral is taken between the bounce points (Northrop & Teller Reference Northrop and Teller1960). The rapid bouncing enforces parity symmetry along the magnetic-field direction. A direct result of both
$\mu$
and
$J$
conservation is then (2.38). The correspondence can be seen more directly by combining (2.27), (2.33), (2.37), and (2.38) to obtain
which is consistent with the conservation of
$J$
as defined in (2.48).
Thus, the conservation of
$\mu$
and
$J$
in magnetised, collisionless plasmas provides a concrete example in which the assumptions made in § 2.2 hold, and (2.41) and (2.42) yield accurate EoSs. However, the general approach we have presented here shows that adiabaticity in the fluid-dynamical context arises due to specific symmetries and phase-space conservation properties of the distribution function. It is best highlighted by our isotropic example in § 2.1 that individual particles are not required to possess adiabatic invariants for adiabatic EoSs to hold, although that is the case for rapidly gyrating, trapped particles in magnetised collisionless plasmas. Consequently, the requirement of particle trapping is not necessary for the double-adiabatic equations (2.41) and (2.42) to be valid. The EoSs would hold as long as our assumptions about the symmetries of
$f$
and conservation of phase-space volume and area are true. This is why, for example, the CGL equations are found to hold in homogeneous, magnetised, collisionless plasmas, where trapped particles do not exist.
3. Double-adiabatic equations of state for ultra-relativistic plasmas
To make analytical progress in evaluating (2.41) and (2.42), we now consider the ultra-relativistic limit,
$\tilde {p}_\perp , \tilde {p}_\parallel \gg mc$
. For ease of notation, we drop the tildes above the rest-frame quantities. We employ the coordinate transformation
$p_\perp = p \sqrt {1 - \xi ^2}$
and
$p_\parallel = p \xi$
, where
$\xi$
is the cosine of the pitch angle with respect to the magnetic field. The pressure-evolution equations (2.41) and (2.42) then simplify to
The ‘total’ pressure (2.21), here equal to
$(2P_\perp + P_\parallel )/3$
, is given by
\begin{align} P = \frac {2c \pi }{3} &\frac {n^{\prime 2}}{B'}\int _{-1}^1 {\rm d}\xi \int _0^\infty {\rm d}p \nonumber \\[4pt] & \times p^3[\xi ^2 + (1-\xi ^2)B^{\prime 3}/n^{\prime 2}]^{1/2} f_0(p, \xi ), \end{align}
and has the value
at
$n' = B' = 1$
.
If the initial distribution function is isotropic, i.e.
$f_0(p, \xi ) = f_0(p)$
, the integrals over the momentum magnitude can be replaced by (3.4). In this limit, the expressions for
$P_\perp$
and
$P_\parallel$
are
where
\begin{align} I_\perp (A) &= -\frac {3}{4} \left [\frac {1}{1 - A} + \frac {(A - 2)}{(1 - A)^{3/2}}\,\operatorname {arsinh}\!\left (\sqrt {\frac {1}{A} - 1}\right )\right ],\\[-10pt]\nonumber \end{align}
\begin{align} I_\parallel (A) &= \frac {3}{2} \left [\frac {1}{1 - A} - \frac {A}{(1 - A)^{3/2}}\,\operatorname {arsinh}\!\left (\sqrt {\frac {1}{A} - 1}\right )\right ], \end{align}
and
Here, we have also defined
$P_\perp ' \equiv P_\perp / P_{0}$
and
$P_\parallel ' \equiv P_\parallel / P_{0}$
. In the case
$A \gt 1$
, the identity
$\operatorname {arsinh}\xi = -i\arcsin (i\xi )$
can be used to rewrite (3.7) and (3.8) in an explicitly real form. Key dimensionless physical parameters, viz.
$\beta _\perp \equiv 8\pi P_\perp /B^2$
,
$\beta _\parallel \equiv 8\pi P_\parallel /B^2$
, and the pressure anisotropy
$P_\perp /P_\parallel$
, depend only on the parameter
$A$
. For example, the latter is
which we plot in figure 1.
Ratio
$P_\perp /P_\parallel$
as a function of
$A \equiv B^{\prime 3}/n^{\prime 2}$
(see (3.10)), shown by the green line. The solid light- and dark-grey lines indicate
$P_\perp /P_\parallel = A$
and
$P_\perp /P_\parallel = A^{4/5}$
, respectively. In the limits
$P_\perp /P_\parallel \ll 1$
and
$P_\perp /P_\parallel \gg 1$
, the ratio scales as
$P_\perp /P_\parallel \sim A$
, whereas close to isotropy it scales as
$P_\perp /P_\parallel \sim A^{4/5}$
.

The relations (3.5) and (3.6) are the double-adiabatic EoSs for an ultra-relativistic collisionless magnetised plasma with an initially isotropic distribution function. As the point
$t=0$
is arbitrary, the above EoS holds as long as
$f$
is isotropic at some point in time. An appealing property of these relations is their independence of the specific form of
$f_0(p)$
, provided that most particles contributing to the pressure tensor are ultra-relativistic. If we impose isotropic evolution,
$P_\perp = P_\parallel \equiv P$
, then setting (3.1) and (3.2) equal to each other yields
$B^{\prime 3} = n^{\prime 2}$
. Inserting this back into either of the pressure equations gives the expected result for an isotropic ultra-relativistic plasma
It is useful to consider the EoSs (3.5) and (3.6) in three key asymptotic limits:
-
(i) large perpendicular pressure,
$P_\perp \gg P_\parallel \implies A\gg 1$
:(3.12)
\begin{align} P_\perp ' = \frac {3\pi }{8}n' B^{\prime 1/2}, \quad P_\parallel ' =\frac {3 \pi }{4} \frac {n^{\prime 3}}{B^{\prime 5/2}}; \end{align}
-
(ii) almost isotropic plasma,
$P_\perp \approx P_\parallel \implies |A - 1| \ll 1$
:(3.13)
\begin{align} P_\perp ' = (n'B')^{4/5}, \quad P_\parallel ' = \left (\frac {n{'}^{3}}{B{'}^{2}}\right )^{4/5}; \end{align}
-
(iii) large parallel pressure,
$P_\perp \ll P_\parallel \implies A\ll 1$
:(3.14)
\begin{align} P_\perp ' = \frac {3}{4} B^{\prime 2}\ln \left ( \frac {n^{\prime 2}}{B^{\prime 3}}\right ), \quad P_\parallel ' = \frac {3}{2}\frac {n^{\prime 2}}{B'}. \end{align}
The asymptotic limits of significant pressure anisotropy agree with the preliminary, heuristic discussion in Gedalin & Oiberman (Reference Gedalin and Oiberman1995), except for the logarithmic correction in (3.14), which arises from a careful treatment of the limit
$P_\perp \ll P_\parallel$
. In such highly anisotropic regimes, the EoSs can be derived straightforwardly from scaling arguments. For example, when
$P_\perp \gg P_\parallel$
, we have
$\gamma \propto p_\perp$
. From the conservation of
$\mu$
(2.47), the perpendicular momentum scales as
$p_\perp \propto B^{1/2}$
. Conservation of
$J$
(2.48) implies
$p_\parallel \propto n/B$
. Substituting these scalings into (2.40), we find
$P_\perp \propto n B^{1/2}$
and
$P_\parallel \propto n^{3}/B^{5/2}$
. A similar argument applies in the limit
$P_\parallel \gg P_\perp$
, where
$\gamma \propto p_\parallel$
, yielding
$P_\perp \propto B^2$
and
$P_\parallel \propto n^2/B$
. This discussion also implies that the EoSs obtained in the asymptotic limits
$P_\perp \gg P_\parallel$
and
$P_\parallel \gg P_\perp$
hold even if
$f_0$
is not isotropic, as found in Gedalin & Oiberman (Reference Gedalin and Oiberman1995), and as can be shown directly from (3.1) and (3.2).
3.1. Generalised adiabatic indices
For ease of use in plasma fluid models, we express
$I_\perp$
and
$I_\parallel$
in the form of power laws,
$I_\perp (A) = A^{\varGamma _\perp {-2}}$
and
$I_\parallel (A) = A^{\varGamma _\parallel {-1}}$
, where
$\varGamma _\perp$
and
$\varGamma _\parallel$
are, in general, functions of
$A$
. The pressure-evolution equations (3.5) and (3.6) are then
We refer to
$\varGamma _\perp$
and
$\varGamma _\parallel$
as the perpendicular and parallel generalised adiabatic indices, respectively, and plot them as a function of
$A$
in figure 2. The values of
$\varGamma _\perp$
and
$\varGamma _\parallel$
in the asymptotic limits of
$A$
can be found by comparing (3.15) with (3.12), (3.13) and (3.14). Close to isotropy,
$|A-1|\ll 1$
,
$(\varGamma _\perp , \varGamma _\parallel ) = ({8}/5, {4}/5)$
and for
$A\gg 1$
,
$(\varGamma _\perp , \varGamma _\parallel ) = ({3}/2, 1/2)$
. When
$A\ll 1$
,
$\varGamma _\parallel = {1}$
, but
$\varGamma _\perp$
tends to 2 only slowly due to the logarithmic dependence of
$I_\perp (A)$
on
$A$
in this limit.
Generalised adiabatic indices
$\varGamma _\perp$
(left) and
$\varGamma _\parallel$
(right), defined in § 3.1, as functions of the pressure anisotropy parameter
$A = B^{\prime 3}/n^{\prime 2}$
. The pressure-evolution equations written in terms of these generalised adiabatic indices are (3.15). The dashed vertical lines correspond to the point
$A= 1$
, i.e. isotropic plasma.

3.2. Eulerian formulation
We have so far worked in a Lagrangian framework with respect to the bulk velocity
$\boldsymbol{u}$
. Numerical implementation of (3.1) and (3.2), however, is often simpler when the EoSs are written in their Eulerian form
Here,
${\rm d}/{\rm d}t \equiv \partial /\partial t + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla}$
is the convective time derivative, and
$A = A(P_\perp /P_\parallel )$
is related implicitly to the pressures by (3.10). The absence of primed quantities makes these equations straightforward to implement in fluid simulations and facilitates the inclusion of additional non-ideal effects, such as cooling.
4. PIC simulations of ultra-relativistic plasma in a compressing box
To test the double-adiabatic EoSs that we have derived, we carry out fully kinetic, two-dimensional PIC simulations using the relativistic code OSIRIS (Fonseca et al. Reference Fonseca2002). We employ a set-up in which an electron–positron plasma is compressed perpendicularly to a straight uniform magnetic field of initial strength
$B_0$
. As a result, both the density and magnetic field increase while maintaining
$B' = n'$
, meaning that, in this set-up,
$A = n'$
. We initialise each species with a Maxwell–Jüttner distribution
with temperature
$T_0 = 10 \, m_e c^2$
and density
$n_0$
such that
$\beta _0 = 8 \pi n_0T_0/B_0^2= 1$
for both species. We choose these parameters in order to postpone the onset of pressure-anisotropy-driven kinetic instabilities (Chandrasekhar et al. Reference Chandrasekhar, Kaufman and Watson1958; Hasegawa Reference Hasegawa1969; Southwood & Kivelson Reference Southwood and Kivelson1993; Bott et al. Reference Bott, Cowley and Schekochihin2024), which will eventually break the double-adiabatic evolution. The full details of our simulations are discussed in Appendix A.
We observe that, as a result of the compression, both the true and modified pressures increase with time as predicted. Plots of these quantities obtained from our simulations, as functions of
$n'$
, are shown in the left panel of figure 3. We find excellent agreement with our theoretical predictions. The evolution of
$P_\perp$
and
$P_\parallel$
satisfies (3.5) and (3.6) with
$B' = n'$
. Figure 3 shows that at values of
$n'$
close to unity, i.e. for small anisotropies,
$P'_\perp = n^{\prime 8/5}$
and
$P'_\parallel = n^{\prime 4/5}$
, as expected from (3.13), until
$n' \approx 3.2$
, when the assumption
$P_\perp \sim P_\parallel$
no longer holds and the full pressure equations are required. Our simulations also show that the evolution of
$\hat {P}'_\perp$
and
$\hat {P}'_\parallel$
is precisely captured by (2.46).
The true (2.40) and modified (2.43) pressures, calculated from PIC simulations of plasma undergoing compression perpendicular (left panel) and parallel (right panel) to the magnetic field, plotted as functions of the normalised density
$n'$
. In the perpendicular-compression case,
$ B' = n'$
, while in the parallel case
$ B' = 1$
. The theoretical predictions given by (3.5), (3.6), and (2.46) are shown as solid grey lines. In the left panel, the grey dashed lines correspond to
$ P_\perp ' = n^{\prime 8/5}$
and
$ P_\parallel ' = n^{\prime 4/5}$
, and in the right panel to
$ P_\perp ' = n^{\prime 4/5}$
and
$ P_{\parallel }' = n^{\prime 12/5}$
. At larger compressions, a clear departure from double-adiabatic theory is observed due to the excitation of the relevant pressure-anisotropy-driven kinetic instabilities.

Figure 3. Long description
Two line graphs depict the true and modified pressures in plasma undergoing compression perpendicular and parallel to the magnetic field. Panel A: The left panel shows the perpendicular compression. The x-axis is labeled ‘log n’’ ranging from 0 to 1.05, and the y-axis is labeled with a range from 0 to 1.5. The graph includes four data series: log P’ with blue circles, log P’ with red circles, log P’ with green circles, and log P’’ with purple circles. Solid grey lines represent theoretical predictions. Grey dashed lines correspond to specific theoretical values. Panel B: The right panel shows the parallel compression. The x-axis is labeled ‘log n’’ ranging from 0 to 0.35, and the y-axis is labeled with a range from 0 to 1. The graph includes the same four data series as the left panel. Solid grey lines represent theoretical predictions. Grey dashed lines correspond to specific theoretical values. Both panels show a clear departure from double-adiabatic theory at larger compressions due to the excitation of relevant pressure-anisotropy-driven kinetic instabilities.
Eventually, when
$\varDelta \sim 1/\beta _\perp$
, around
$n' \approx {6.3}$
, we begin to observe an increase in magnetic energy, indicating the onset of the relativistic mirror instability (Galishnikova, Philippov & Quataert Reference Galishnikova, Philippov and Quataert2023). As a result of the microscale fluctuations in the magnetic field excited by the instability, double-adiabatic evolution breaks down at
$n' \approx {7.9}$
due to
$\mu$
no longer being conserved. Hence, the pressures and modified pressures cease to satisfy (3.5), (3.6), and (2.46). In future work, we will further study the late-time evolution of these systems, where the onset of kinetic instabilities becomes important.
To check the general form of the ultra-relativistic pressure-evolution equations (3.5) and (3.6), we also ran simulations in which the plasma was compressed parallel to the magnetic field. The simulation results are shown in the right panel of figure 3. During such compression,
$n'$
still increases with time, but
$B'$
remains fixed at unity. We again found good agreement with our double-adiabatic evolution equations, observing initial growth of the pressures given by
$P_\perp ' = n^{\prime 4/5}$
and
$P_\parallel ' = n^{\prime 12/5}$
, as predicted by (3.13). In this arrangement,
$A = 1/n^{\prime 2}$
, and so during the compression a negative pressure anisotropy arises. As a result, we observe the onset of the firehose instability at
$\varDelta \sim -1/\beta _\parallel$
, occurring at
$n' \approx 1.9$
. Once the firehose fluctuations become sufficiently large, they begin to scatter particles, thereby double-adiabatic evolution no longer holds due to a breaking of
$\mu$
conservation. This behaviour is evident in figure 3 at
$n' \approx 2.1$
, where the pressure and modified-pressure evolutions depart from (3.5), (3.6), and (2.46).
5. Discussion
In this paper, we have shown that phase-space symmetries naturally lead to adiabatic evolution equations for the pressure(s). We did so by treating the distribution function as the density of a six-dimensional, incompressible fluid and demonstrating that imposing such symmetries leads to a self-similar evolution of the distribution function. Taking the appropriate moments of the evolved distribution function then yields adiabatic evolution equations whose form, in general, depends on the initial distribution function.
For a distribution function that is isotropic in the peculiar momentum, such as that of a plasma with efficient pitch-angle scattering, we found an adiabatic EoS for the pressure, which reduces to
$P \propto n^{5/3}$
and
$P \propto n^{4/3}$
in the non-relativistic and ultra-relativistic limits, respectively. We also considered a distribution whose peculiar-momentum dependence is gyrotropic with respect to a vector field,
$\boldsymbol{B}$
(the magnetic field), frozen into the bulk flow and symmetric under parity along that direction. We found that such a system satisfies double-adiabatic equations. A natural realisation of this scenario is a collisionless, magnetised plasma, in which particles gyrate rapidly around the magnetic field and undergo fast bounce motion in magnetic mirrors. In this case, the double-adiabatic evolution equations govern the parallel and perpendicular pressures, and our symmetry assumptions are consistent with the conservation of the first and second adiabatic invariants,
$\mu$
and
$J$
.
Through developing this symmetry-based formalism, we extended the double-adiabatic evolution equations of a collisionless plasma to the relativistic regime. For an initially isotropic, ultra-relativistically hot plasma, we found
$P_\perp \propto B^2 I_\perp$
and
$P_\parallel \propto n^2 I_\parallel /B$
, where
$I_\perp$
and
$I_\parallel$
are functions of the compression ratio
$B^{\prime 3}/n^{\prime 2}$
, or, equivalently, of the pressure anisotropy. We evaluated these functions in three asymptotic limits and derived the corresponding pressure-evolution equations. Close to isotropy,
$P_\perp \propto (nB)^{4/5}$
and
$P_\parallel \propto (n^{3}/B^{2})^{4/5}$
. When
$P_\perp \gg P_\parallel$
, we obtained
$P_\perp \propto n B^{1/2}$
and
$P_\parallel \propto n^3/B^{5/2}$
. In the opposite limit,
$P_\parallel \gg P_\perp$
, we found
$P_\perp \propto B^2 \ln (n^{\prime 2}/B^{\prime 3})$
and
$P_\parallel \propto {n^2/B}$
. We also expressed our EoSs in a form ready to be implemented in numerical fluid simulations. As the ‘initial’ time is arbitrary, these EoSs hold if the plasma was, at some point in time, isotropic. We tested these pressure-evolution equations numerically using compressing-box PIC simulations and found excellent agreement with the theoretical predictions.
A limitation of the double-adiabatic EoSs presented here is that, for collisionless plasmas, the distribution function is only gyrotropic if the coordinate
$\boldsymbol{r}$
refers to the gyrocentre of the particle. In this case, the statement that the averaged particle velocity is equal to the bulk velocity,
$\langle \dot {\boldsymbol{r}}\rangle _{{t_{\text{fast}}}} = \boldsymbol{u}$
, neglects magnetic drifts and is therefore valid only in the limit of vanishing gyroradius. Additionally, double adiabaticity requires neglect of heat fluxes and collisions. Incorporating such non-ideal and finite-gyroradius effects would therefore require corrections to the pressure-evolution equations derived here.
Our derivation was also restricted to plasmas with a non-relativistic bulk flow. An extension to the fully relativistic regime is possible and yields the same EoSs, with the time coordinate replaced by the proper time of the bulk flow (Wierzchucka et al. Reference Wierzchucka, Bilbao, Uzdensky and Schekochihin2026b
). In this regime,
$n$
and
$B$
denote, respectively, the density and magnetic-field strength in the comoving frame. Constructing a covariant Lagrangian description of phase-space evolution is mathematically involved, owing, e.g., to the need for relativistic Clebsch variables. It turns out to be considerably easier to derive these results from the relativistic drift-kinetic equation (Wierzchucka et al. Reference Wierzchucka, Bilbao, Uzdensky and Schekochihin2026b
).
The physical set-up that we consider here provides a natural starting point for a study of driven kinetic processes, such as the excitation of the firehose and mirror instabilities, in high-energy plasmas. Existing studies suggest that relativistic effects introduce order-unity modifications to the instability thresholds (Galishnikova et al. Reference Galishnikova, Philippov and Quataert2023; Zhdankin et al. Reference Zhdankin, Kunz and Uzdensky2023), yet many aspects of their nonlinear evolution and saturation in the relativistic regime remain poorly understood, in particular, how driven relativistic mirror and firehose instabilities evolve in the nonlinear regime. The relativistic EoSs derived in this work already indicate departures from non-relativistic behaviour. For example, the perpendicular-pressure evolution (3.5) implies that, unlike in the non-relativistic regime,
$\beta _\perp = 8\pi P_\perp / B^2$
decreases under compression perpendicular to the magnetic field. As a result, significantly larger compressions are required to trigger instabilities associated with positive pressure anisotropy, suggesting that relativistic collisionless plasmas could be more stable to anisotropy-driven instabilities than their non-relativistic counterparts.
Our theory of double-adiabatic evolution is relevant to other processes occurring in relativistic plasmas, such as the emission of synchrotron radiation. Recent studies have found that synchrotron radiation generates a negative pressure anisotropy in a cooling plasma, exciting the firehose and maser instabilities (Bilbao & Silva Reference Bilbao and Silva2023; Zhdankin et al. Reference Zhdankin, Kunz and Uzdensky2023; Bilbao et al. Reference Bilbao, Ewart, Assunçao, Silva and Silva2024). When the firehose instability interacts with the synchrotron cooling instability (Simon & Axford Reference Simon and Axford1967), a two-phase medium forms (Wierzchucka et al. Reference Wierzchucka, Bilbao, Ewart, Uzdensky and Schekochihin2026a
). One of the phases is high-
$\beta$
, and its dynamics is dominated by the firehose instability, which keeps the pressure anisotropy pinned to the firehose marginality condition. The second phase is low-
$\beta$
, and therefore no firehose fluctuations exist within it. The EoSs found in this work have enabled an effective description of such a medium using a fluid model in which the two phases have different effective adiabatic indices (Wierzchucka et al. Reference Wierzchucka, Bilbao, Ewart, Uzdensky and Schekochihin2026a
). The determination of effective adiabatic indices through double-adiabatic EoSs is also important for jump conditions of magnetised, collisionless shocks (Sironi & Spitkovsky Reference Sironi and Spitkovsky2009; Sironi, Spitkovsky & Arons Reference Sironi, Spitkovsky and Arons2013; Bret & Narayan Reference Bret and Narayan2018).
This work may be important for the study of plasmoid dynamics in magnetic reconnection in both astrophysical (Øieroset et al. Reference Øieroset, Phan, Fujimoto, Lin and Lepping2001; Retinò et al. Reference Retinò, Sundkvist, Vaivads, Mozer, André and Owen2007) and laboratory (Hare et al. Reference Hare2017; Totorica, Abel & Fiuza Reference Totorica, Abel and Fiuza2017) settings. Reversing magnetic-field configurations are unstable and collapse to form current sheets (Furth, Killeen & Rosenbluth Reference Furth, Killeen and Rosenbluth1963). At sufficiently high Lundquist number, these sheets fragment into secondary magnetic islands, or plasmoids (Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010). In the relativistic regime, magnetic reconnection has been observed to generate extended non-thermal distributions and accelerate particles (Sironi & Spitkovsky Reference Sironi and Spitkovsky2014; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2015). Several mechanisms have been proposed to explain this acceleration and the formation of power-law tails, including plasmoid mergers and Fermi-type processes (Drake et al. Reference Drake, Swisdak, Che and Shay2006; Li et al. Reference Li, Guo, Li and Li2018). Double-adiabatic models based on the conservation of
$\mu$
and
$J$
have been used to describe Fermi acceleration, although primarily in the non-relativistic regime (Le et al. Reference Le, Egedal, Daughton, Fox and Katz2009; Montag et al. Reference Montag, Egedal, Lichko and Wetherton2017). PIC simulations of relativistic reconnection have also shown that conservation of the first two adiabatic invariants can energise particles within compressing plasmoids during the ‘secondary acceleration phase’, following initial acceleration in current sheets (Petropoulou & Sironi Reference Petropoulou and Sironi2018; Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021). The double-adiabatic EoSs derived here offer a useful model for describing particle energisation within plasmoids formed during relativistic reconnection.
To conclude, in this paper, we have used a symmetry-based approach to extend the double-adiabatic equations to the relativistic regime, providing a ready-to-implement fluid description for collisionless, magnetised relativistic plasmas. The generality of the formalism suggests that it may also be applicable to other systems satisfying similar assumptions, such as the presence of a frozen-in vector field.
Note: right before the submission of this manuscript, we became aware of similar results that were obtained concurrently and independently by Ley, Tran & Zweibel (Reference Ley, Tran and Zweibel2026) and submitted simultaneously with this paper. There are two key differences between our work and that of Ley et al. (Reference Ley, Tran and Zweibel2026). First, they focus on deriving the adiabatic evolution of
$f$
using the (homogeneous) drift-kinetic equation, as opposed to phase-space symmetries. Secondly, while we show that the relativistic double-adiabatic equations (3.5) and (3.6) hold for any distribution function that was at some point in time isotropic, Ley et al. (Reference Ley, Tran and Zweibel2026) limit themselves to an initial Maxwell–Jüttner equilibrium. Additionally, also shortly before submission, we learned about an otherwise unpublished PhD thesis by P. Simeon (Reference Simeon2014) containing a preliminary calculation of adiabatic indices in relativistic, collisionless plasmas. In one chapter, the expressions (2.41) and (2.42) are derived, and the integrals are evaluated in the limit close to the initial state, i.e.
$B' \approx 1$
and
$n' \approx 1$
, yielding the EoSs (3.13).
Acknowledgements
We thank A.F.A. Bott, P.G. Ivanov, E.J. Kolmes, M. Lyutikov and M.L. Nastac for useful discussions. We are especially thankful to F. Bacchini for giving us access to the compressing-box module in the Zeltron code, which inspired us to pursue this project.
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
Simulations were performed at Deucalion (Portugal), funded by FCT MAPs 2: Masers in Astrophysical Plasmas I.P. project 2025.00058.CPCA.A3, and at MareNostrum 5 (Spain), funded by EuroHPC-JU extreme scale access project KEEP: Kinetic Electrodynamics in Extreme Plasmas (EHPC-EXT-2025E02-066). A.W. was supported by a Clarendon Scholarship and the Merton College Tira Wannamethee Scholarship. P.J.B. was supported by a Leverhulme–Peierls Fellowship. A.G.R.T. was supported by US DOE grant DE-SC0024562 and NNSA Center of Excellence DE-NA0003869. The work of A.A.S. and P.J.B. was supported in part by the UK STFC grant ST/W000903/1. A.A.S. and A.W. were also supported in part by the Simons Foundation via the Simons Investigator Award to A.A.S. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Details of PIC simulations
To investigate numerically the double-adiabatic EoSs of relativistic plasmas, we use OSIRIS (Fonseca et al. Reference Fonseca2002), a fully relativistic, three-dimensional, highly parallelised PIC code. The PIC numerical method works by solving the equations of motion for charged particles subject to the Lorentz force and by calculating the self-consistent fields on a grid. It is commonly used in plasma physics because of its ability to model particle interactions accurately, thereby allowing the capture of effects beyond the fluid approximation.
We initialised our simulations with a magnetic field
$\boldsymbol{B}_0 = B_0\hat {\boldsymbol{x}}$
and allowed spatial variations only in the (
$x$
,
$y$
) plane. This choice reflects the fact that kinetic processes, such as the mirror instability, produce fluctuations with wave vectors both parallel and perpendicular to the magnetic field. Our plasma consisted of electrons and positrons of mass
$m_e$
. At the start of the simulation, each species had a Maxwell–Jüttner distribution (4.1) with temperature
$T_0 = 10 \, m_e c^2$
and an associated (ultra-relativistic) gyration frequency
$\varOmega _{e0} = e B_0/ 3 \theta _0 m_e {c}$
, where
$\theta _0 \equiv T_0/m_ec^2$
. All times and distances were normalised to the initial gyration period
$\varOmega _{e0}^{-1}$
and the Larmor radius
$\rho _{e0} = c/\varOmega _{e0}$
, respectively. We set the initial densities of the positrons and electrons to
$n_0 = \omega _{pe}^2 m_e/ 4\pi e^2$
, where
$\omega _{pe}$
is the (non-relativistic) plasma frequency. As we are interested in purely collisionless dynamics undisturbed by instabilities, we used
$\beta _0 = 8\pi n_0 T_0/B_0^2 = 1$
, corresponding to
$B_0 = 4.47 \ \omega _{pe} m_e c/e$
and
$\omega _{pe}/\varOmega _{e0} = {6.7}$
.
Since the EoSs (3.1) and (3.2) express perpendicular and parallel pressures in terms of the density and magnetic-field strength, we set up a large-scale compressing flow in our simulations to introduce controlled variations in
$B'$
and
$n'$
. The standard PIC treatment of a background, non-relativistic compressing flow, described in Sironi & Narayan (Reference Sironi and Narayan2015), is to solve the PIC equations in a frame comoving with the compressing flow. We implemented this method in OSIRIS by modifying the existing shearing-box algorithm (Inchingolo et al. Reference Inchingolo, Grismayer, Loureiro, Fonseca and Silva2018). For a general compression, the laboratory coordinate
$\boldsymbol{x}$
is related to the comoving coordinate
$\boldsymbol{x'}$
by
$\boldsymbol{x} = \boldsymbol{L} \boldsymbol{x'}$
, where
\begin{align} \boldsymbol{L} = \frac {\partial \boldsymbol{x}}{\partial \boldsymbol{x'}} = \begin{pmatrix} q_x & \quad 0 & \quad 0 \\[4pt] 0 & \quad q_y & \quad 0 \\[4pt] 0 & \quad 0 & \quad q_z \end{pmatrix}, \end{align}
and
$q_x$
,
$q_y$
and
$q_z$
are time-dependent compression coefficients. To model a compression perpendicular to the magnetic field, we set
where
$t$
is time and
$S$
is the compression rate. In order for the magnetic field to vary slowly compared with the gyration time, we choose
$S/\varOmega _{e0} = 10^{-2}$
. As a result of this compression, the magnetic field and density satisfy
$B' = n' = (1 + St)^2$
. Our set-up is summarised in figure 4. We also ran simulations in which the plasma was compressed parallel to the guide field, i.e.
with all other parameters identical to those of the perpendicular-compression simulation. In this case,
$n' = 1 + St$
, whereas the magnetic field remains constant during the compression.
Our simulation set-up. The background magnetic field is
$\boldsymbol{B}_0 = B_0 \hat {\boldsymbol{x}}$
, the box is compressed in the
$z$
and
$y$
directions. Only variations in the
$(x,y)$
plane are allowed in the simulation.

Our simulations used a spatial resolution
$\Delta x = 0.0{36} \, \rho _{e0}$
in each direction and a time step of
$\Delta t = 0.00{55} \, \varOmega _{e0}^{-1}$
. This ensured that the Courant–Friedrichs–Lewy condition in two dimensions,
$c \Delta t / \Delta x \lt {1/}\sqrt {2}$
, was satisfied throughout the entire duration of the compression and that Larmor-scale motions were properly resolved. The domain size was
$L_x \times L_y = 20{.1} \rho _{e0} \times 20.{1} \rho _{e0}$
, with a corresponding number of cells
$N_x \times N_y = 552 \times 552$
. We imposed periodic boundary conditions in both directions and used
$34{2}$
particles per cell with cubic particle shape interpolation.



P⊥/P∥
A≡B′3/n′2
P⊥/P∥=A
P⊥/P∥=A4/5
P⊥/P∥≪1
P⊥/P∥≫1
P⊥/P∥∼A
P⊥/P∥∼A4/5
Γ⊥
Γ∥
A=B′3/n′2
A=1
n′
B′=n′
B′=1
P⊥′=n′8/5
P∥′=n′4/5
P⊥′=n′4/5
P∥′=n′12/5
B0=B0x^
z
y
(x,y)