Synchrotron Firehose Instability

We consider a population of electrons (or positrons) with an ultra-relativistic temperature (such that momenta p ≫ m_e c) being cooled by synchrotron radiation. The classical synchrotron radiation reaction force can be expressed as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{sync}}=-\displaystyle \frac{2}{3}{r}_{e}^{2}{\gamma }^{2}\left[{\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}}{c}\times {\boldsymbol{B}}\right)}^{2}-{\left(\displaystyle \frac{{\boldsymbol{v}}\cdot {\boldsymbol{E}}}{c}\right)}^{2}\right]\displaystyle \frac{{\boldsymbol{v}}}{c},\end{eqnarray} \tag{ 1 }$$

where r_e ≡ e²/m_e c² is the classical electron radius, v = p /(γ m_e ) is the particle velocity, and $\gamma \equiv {\left(1-{v}^{2}/{c}^{2}\right)}^{-1/2}$ is the Lorentz factor (Landau & Lifshitz 1975; Rybicki & Lightman 1991). Specializing to the case of a uniform magnetic field of B = B ₀ and a vanishing electric field of E = 0, Equation (1) can be conveniently written as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{sync}}=-\displaystyle \frac{2}{3}\displaystyle \frac{{r}_{e}^{2}{B}_{0}^{2}}{{m}_{e}^{2}{c}^{2}}{p}^{2}{\sin }^{2}\theta \,\hat{{\boldsymbol{p}}}=-\,\displaystyle \frac{{p}^{2}{\sin }^{2}\theta }{{\overline{p}}_{0}{\tau }_{\mathrm{cool}}}\,\hat{{\boldsymbol{p}}},\end{eqnarray} \tag{ 2 }$$

where $\theta \equiv {\cos }^{-1}({\boldsymbol{B}}\cdot {\boldsymbol{p}}/{Bp})$ is the particle's pitch angle, ${\overline{p}}_{0}$ is the initial average momentum, and

$$\begin{eqnarray}&&{\tau }_{\mathrm{cool}}\equiv \displaystyle \frac{3{m}_{e}^{2}{c}^{2}}{2{r}_{e}^{2}{B}_{0}^{2}{\overline{p}}_{0}}\end{eqnarray} \tag{ 3 }$$

is the characteristic cooling timescale for the average particle.

The evolution of a uniform, gyrotropic electron distribution f_e ( p , t) = f_e (p, θ, t) subject to the force in Equation (2) is described by the Vlasov equation

$$\begin{eqnarray}&&{\partial }_{t}{f}_{e}=-\displaystyle \frac{\partial }{\partial {\boldsymbol{p}}}\,\cdot \,\left({{\boldsymbol{F}}}_{\mathrm{sync}}{f}_{e}\right)=\displaystyle \frac{{\sin }^{2}\theta }{{\tau }_{\mathrm{cool}}{\overline{p}}_{0}{p}^{2}}\,{\partial }_{p}\left({p}^{4}{f}_{e}\right).\end{eqnarray} \tag{ 4 }$$

Given an isotropic initial distribution f_e0(p), this equation can be readily solved by the method of characteristics to obtain

$$\begin{eqnarray}&&{f}_{e}(p,\theta ,t)=\displaystyle \frac{{f}_{e0}\left(p/[1-{pt}{\sin }^{2}\theta /({\overline{p}}_{0}{\tau }_{\mathrm{cool}})]\right)}{{\left[1-{pt}{\sin }^{2}\theta /({\overline{p}}_{0}{\tau }_{\mathrm{cool}})\right]}^{4}}\end{eqnarray} \tag{ 5 }$$

for $p\lt {p}_{\mathrm{cut}}(\theta )\equiv {\overline{p}}_{0}{\tau }_{\mathrm{cool}}/(t{\sin }^{2}\theta )$ and f_e = 0 for p ≥ p_cut(θ). For an initial ultra-relativistic Maxwell–Jüttner distribution, ${f}_{e0}={n}_{0}\exp (-p/{p}_{T})/(8\pi {p}_{T}^{3})$, Equation (5) becomes

$$\begin{eqnarray}&&{f}_{e}(p,\theta ,t)=\displaystyle \frac{{n}_{0}}{8\pi {p}_{T}^{3}}\displaystyle \frac{\exp \left[\tfrac{-p/{p}_{T}}{1-{pt}{\sin }^{2}\theta /(3{p}_{T}{\tau }_{\mathrm{cool}})}\right]}{{\left[1-{pt}{\sin }^{2}\theta /(3{p}_{T}{\tau }_{\mathrm{cool}})\right]}^{4}},\end{eqnarray} \tag{ 6 }$$

where ${p}_{T}\equiv {T}_{e0}/c={\overline{p}}_{0}/3$ is the thermal momentum and n₀ is the number density. The solutions given by Equations (5)–(6) have a hard cutoff at p = p_cut(θ), due to high-energy particles catching up to low-energy particles as they cool (causing a "phase-space shock"). The maximum of the distribution in Equation (6) at a given momentum can be shown (from ∂f_e /∂θ = 0) to occur at a pitch angle satisfying ${\sin }^{2}\theta =3({p}_{T}/p-1/4)({\tau }_{\mathrm{cool}}/t)$ for ${(1/4+t/3{\tau }_{\mathrm{cool}})}^{-1}\lt p/{p}_{T}\lt 4$, at θ = 0 for p/p_T > 4, and at θ = π/2 for $p/{p}_{T}\lt {(1/4+t/3{\tau }_{\mathrm{cool}})}^{-1};$ thus high-energy particles are focused in the direction of B ₀. Iso-contours of the solution given by Equation (6) at t = τ_cool are plotted in the top panel of Figure 1. In the middle panel of Figure 1, we show the distribution averaged over angles, ${p}^{2}{f}_{e}(p,t)\equiv (1/2){p}^{2}\int d\cos \theta \,{f}_{e}(p,\theta ,t)$, at times t/τ_cool ∈ {0, 0.25, 1, 4, 16}. At late times, the solution develops a power-law range with a scaling of p² f_e (p) ∝ p⁻² from p ∼ p_cut(π/2) to $p\sim {\overline{p}}_{0}$, which is a consequence of the superposition of phase-space shocks at different p_cut(θ). In the bottom panel of Figure 1, we show the evolution of the distribution in the perpendicular direction, p² f_e (p, θ = π/2), highlighting the gradual steepening of the phase-space shock.

Zoom In Zoom Out Reset image size

For early times ≡ t/τ_cool ≪ 1, Equation (6) can be Taylor-expanded as long as the momenta considered are not too large, p ≲ p_T . To first order in , the result is

$$\begin{eqnarray}&&{f}_{e}\approx \displaystyle \frac{{n}_{0}}{8\pi {p}_{T}^{3}}\left[1+\left(4-\displaystyle \frac{p}{{p}_{T}}\right)\displaystyle \frac{{pt}}{3{p}_{T}{\tau }_{\mathrm{cool}}}{\sin }^{2}\theta \right]{e}^{-p/{p}_{T}}.\end{eqnarray} \tag{ 7 }$$

The distribution retains its shape in the parallel direction (θ = 0), while in the perpendicular direction (θ = π/2) it has a pileup of particles for p < 4p_T and a depletion of particles for p > 4p_T .

The evolution of the electron kinetic energy density n₀ E_kin,e and the parallel and perpendicular components of the electron pressure tensor may be obtained by taking the appropriate moments of Equation (7). The results are

$$\begin{eqnarray}\begin{array}{rcl}{n}_{0}{E}_{{\rm{kin}},e} & = & \displaystyle \int {d}^{3}p\,{{pcf}}_{e}\approx 3{n}_{0}{p}_{T}c\left(1-\displaystyle \frac{8}{9}\displaystyle \frac{t}{{\tau }_{{\rm{cool}}}}\right),\\ {P}_{\parallel ,e} & = & \displaystyle \int {d}^{3}p\,{p}_{\parallel }{v}_{\parallel }{f}_{e}\approx {n}_{0}{p}_{T}c\left(1-\displaystyle \frac{8}{15}\displaystyle \frac{t}{{\tau }_{{\rm{cool}}}}\right),\\ {P}_{\perp ,e} & = & \displaystyle \frac{1}{2}\displaystyle \int {d}^{3}p\,{p}_{\perp }{v}_{\perp }{f}_{e}\approx {n}_{0}{p}_{T}c\left(1-\displaystyle \frac{16}{15}\displaystyle \frac{t}{{\tau }_{{\rm{cool}}}}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where the "∥" and "⊥" subscripts indicate directions parallel and perpendicular to B ₀. As a consequence, to first order in , the pressure anisotropy is P_⊥,e /P_∥,e ≈ 1 − (8/15)(t/τ_cool).

In Section 3, we reproduce numerically the analytical solution of Equation (6) by performing a PIC simulation of an uncharged, synchrotron-cooling gas that remains stable to electromagnetic instabilities.

The solution to Equation (6) will be valid until the onset of the firehose instability, which occurs once a sufficient pressure anisotropy is produced to offset the magnetic tension. We anticipate that the instability will grow initially at a super-exponential rate, since the firehose growth rate ${\gamma }_{{\rm{f}}}\equiv d\mathrm{ln}{B}_{{\rm{f}}}/{dt}$ will increase gradually from γ_f = 0 at marginal stability to a value ${\gamma }_{{\rm{f}}}\gtrsim {\tau }_{\mathrm{cool}}^{-1}$ as the plasma evolves into the unstable region of parameter space via continuous cooling and corresponding growth of the pressure anisotropy (unregulated by the instability). Here, we define B_f(t) as the amplitude of the most unstable firehose mode. This super-exponential (or "unregulated") stage may be followed by an exponential (or "regulated") stage once the linear growth rate γ_f becomes large enough that the firehose fluctuations deplete the pressure anisotropy faster than it is replenished by the synchrotron cooling. In practice, the regulated stage may be relatively brief compared to the unregulated stage.

In the Appendix, we calculate analytically the parallel firehose instability criterion from the linearized relativistic Vlasov–Maxwell equations in the limit of long wavelength and low frequency (relative to kinetic scales), using the early-time (t/τ_cool ≪ 1) distribution from Equation (7) as a (slowly varying) background. This instability criterion is

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\perp ,e}}{{P}_{\parallel ,e}}\lt 1-\displaystyle \frac{{C}_{\mathrm{thr}}}{{\beta }_{\parallel ,e}},\end{eqnarray} \tag{ 9 }$$

where C_thr is an order-unity constant and β_∥,e ≡ 8π P_∥,e /B² is the plasma beta using the parallel component of the electron pressure. In our calculation, C_thr = 1 for a relativistic pair plasma (where positrons have the same background cooling distribution as electrons) and C_thr = 2 for an electron–ion plasma (where ions do not cool and thus retain their initial Maxwell–Jüttner distribution f_i0 as the background). Equation (9) agrees with the standard criterion for the relativistic fluid firehose instability (e.g., Barnes & Scargle 1973). The inequality in Equation (9) is satisfied for times after an "onset" time

$$\begin{eqnarray}&&t\gtrsim {\tau }_{\mathrm{onset}}\equiv \displaystyle \frac{15{C}_{\mathrm{thr}}}{8{\beta }_{e0}}{\tau }_{\mathrm{cool}},\end{eqnarray} \tag{ 10 }$$

based on the early-time expressions in Equation (8); this is valid for β_e0 ≫ 1, where β_e0 is the initial electron plasma beta. In terms of physical parameters, the onset time can be expressed as

$$\begin{eqnarray}&&{\tau }_{\mathrm{onset}}=\displaystyle \frac{5{C}_{\mathrm{thr}}}{16}{\left({\theta }_{e0}^{2}{\sigma }_{{\rm{T}}}{n}_{0}c\right)}^{-1},\end{eqnarray} \tag{ 11 }$$

where θ_e0 = T_e0/m_e c² is the dimensionless temperature and ${\sigma }_{{\rm{T}}}=(8\pi /3){r}_{e}^{2}$ is the Thomson cross section. Notably, τ_onset is independent of the magnetic-field strength B₀. For comparison, the (relativistic) Coulomb collision time also scales in proportion to ${({\sigma }_{{\rm{T}}}{n}_{0}c)}^{-1}$, but with a factor of ${\theta }_{e0}^{2}$ rather than ${\theta }_{e0}^{-2}$ (Stepney 1983), so that the onset time is a factor of ${\theta }_{e0}^{-4}$ shorter than the collision time. Thus, we expect Coulomb collisions to be negligible over the onset timescale, as long as the electrons are sufficiently relativistic.

Based on previous studies of the firehose instability in the nonrelativistic regime, we anticipate that finite-Larmor-radius corrections will reduce C_thr by a factor ≈0.7, making the firehose easier to trigger (e.g., Hellinger & Matsumoto 2000; Bott et al. 2021). We also anticipate that, at least in some physical regimes, the firehose instability may occur at oblique rather than parallel orientation (e.g., Li & Habbal 2000).

To illustrate the relevant region of the firehose instability in parameter space, we employ plots of P_⊥,e /P_∥,e versus β_∥,e . In Figure 2, we show the trajectory of the synchrotron-cooling distribution (obtained by integrating Equation (6) numerically) across this parameter space, for initial values of β_∥,e,0 ∈ {2.5, 5, 10, 20, 40} at isotropy. The displayed solutions cover a duration of t/τ_cool ≈ 16, at which time P_⊥,e /P_∥,e ≈ 0.4. We also show the firehose instability thresholds from Equation (9) with C_thr ∈ {0.7, 1, 1.4, 2}. The β_∥,e,0 = 2.5 case misses all of the instability thresholds, while β_∥,e,0 > 10 is sufficient to hit all of the thresholds (after which the development of the SFHI will invalidate the analytical solution). This motivates choosing β_∥,0 = 20 as the fiducial value in the PIC simulations of Section 3, which allows a robust incursion into the unstable zone while avoiding extreme values of β (which are unrealistic and difficult to simulate with PIC).

Zoom In Zoom Out Reset image size

Once the SFHI growth rate becomes positive (at t = τ_onset), we anticipate that the (fast) growth of the most unstable firehose mode at any instant can be modeled using linear theory, but that this growth rate will be modulated on the (slow) cooling timescale due to the continual growth of pressure anisotropy (until regulation). The linear calculation in the Appendix suggests a growth rate of the form

$$\begin{eqnarray}&&{\gamma }_{{\rm{f}}}\equiv \displaystyle \frac{d\mathrm{ln}{B}_{{\rm{f}}}}{{dt}}=\zeta {kc}{\left(\displaystyle \frac{\delta t}{{\tau }_{\mathrm{cool}}}\right)}^{1/2},\end{eqnarray} \tag{ 12 }$$

where ζ is a numerical coefficient that depends on the plasma's composition and temperature (see Equations (A19), (A23), and (A26) in the Appendix) and δ t = t − τ_onset is the time after onset. We can integrate Equation (12) in time to obtain a super-exponential growth

$$\begin{eqnarray}&&{B}_{{\rm{f}}}={B}_{\mathrm{seed}}\exp \left[{\left(\displaystyle \frac{\delta t}{{\tau }_{\mathrm{gr}}}\right)}^{3/2}\right],\end{eqnarray} \tag{ 13 }$$

where ${\tau }_{\mathrm{gr}}={(2\zeta /3)}^{-2/3}{\tau }_{\mathrm{cool}}^{1/3}{({kc})}^{-2/3}$ and B_seed is the initial seed perturbation magnitude (see, e.g., Zhou et al. 2022, for similar arguments applied to the spontaneous development of the Weibel instability). Based on the literature, we expect the most unstable mode to be near the electron kinetic scales; specifically, $k\sim {\rho }_{e0}^{-1}$ for the electron-firehose instability (e.g., Li & Habbal 2000), where ρ_e0 = 3T_e0/eB₀ is the characteristic ultra-relativistic electron Larmor radius. Therefore, ignoring coefficients of order unity, the unregulated growth occurs on a hybrid timescale ${\tau }_{\mathrm{gr}}\sim {\tau }_{\mathrm{cool}}^{1/3}{\tau }_{\mathrm{gyro}}^{2/3}$, where τ_gyro = ρ_e0/c is the electron Larmor timescale. Because the unregulated growth is super-exponential, it is also on this timescale that the pressure anisotropy will be regulated close to the firehose threshold.

To demonstrate the existence of the SFHI and study its properties, we perform a set of radiative PIC simulations using the code Zeltron (Cerutti et al. 2013). The simulations include the synchrotron radiation reaction force, Equation (1), in the particle pusher for electrons and positrons.

Because two spatial dimensions are sufficient to capture the basic features of the SFHI (with either parallel or oblique polarization), we limit ourselves to 2D simulations (an x–y plane with a mean magnetic field of ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{y}}}$ oriented in the domain), which allows us to run simulations with sufficiently long duration to study the slow-cooling regime (τ_gyro/τ_cool ≪ 1). The spatial domain has periodic boundary conditions.

We initialize the PIC simulations with a uniform, isotropic, thermal plasma of electrons and ions (or positrons). In all cases, we initialize both species with the same temperature, T_i0/T_e0 = 1 (where the subscript i refers to either ions or positrons). To study the effect of the plasma's composition on the instability, we consider four cases: (1) an uncharged gas of ultra-relativistic particles, which undergo synchrotron cooling but no electromagnetic interactions, so the instability is artificially inhibited and thus the distribution reproduces the stable analytical solution described in Section 2.1; (2) an ultra-relativistic pair plasma (m_i /m_e = 1) with a dimensionless temperature of θ_e0 ≡ T_e0/m_e c² = 100, where both species radiate; (3) an electron–ion plasma (with "ions" meaning protons, m_i /m_e = 1836) with sub-relativistic ions, taking the dimensionless ion temperature θ_i0 ≡ T_i0/m_i c² = 1/64 ≪ 1, such that electrons are in the ultra-relativistic regime with θ_e0 = (m_i /m_e )θ_i0 ≈ 29 ≫ 1 (we adopt the terminology of Werner et al. 2018 to call this the "semi-relativistic" regime); and (4) an ultra-relativistic electron–ion plasma (m_i /m_e = 1836) with θ_i0 = 100 (and thus θ_e0 ≫ 1 as well).

Our fiducial parameters are ${\beta }_{0}\equiv 8\pi {n}_{0}({T}_{i0}+{T}_{e0})/{B}_{0}^{2}=40$ (such that τ_onset/τ_cool = 3C_thr/32 from Equation (9)) and τ_gyro/τ_cool = 1.2 × 10⁻⁴. The fiducial domain size is L/2π ρ_e0 ≈ 4.1, sufficient to contain the most unstable electron-firehose mode; however, we also performed shorter-duration simulations with L/2π ρ_e0 ≈ 8.1 (which we call "medium size") that give similar results and are used for some of the analysis. We also ran the nonrelativistic ion case on a domain that is four times larger (L/2π ρ_e0 ≈ 16.3 and L/2π ρ_i0 ≈ 3.5, where the nonrelativistic ion gyroradius is ${\rho }_{i0}\approx \sqrt{3{m}_{i}{T}_{i0}}c/{{eB}}_{0}\approx 4.7{\rho }_{e0}$ in this case) with shorter duration, to confirm that ion-scale dynamics do not influence the results. The fiducial simulation lattice is 512² cells (with cell size Δx = ρ_e0/20) and the number of particles per cell per species is 512; the larger cases have the same resolution but correspondingly more cells. The simulations have a timestep of ${\rm{\Delta }}t\approx {\rm{\Delta }}x/(\sqrt{3}c)$ and a fiducial duration of t/τ_cool ∼ 3. A list of the simulations and their primary parameters is shown in Table 1.

In all cases, the plasma initially undergoes a stable cooling phase, in which both the parallel and perpendicular components of the electron (and/or positron) pressure decline, with the perpendicular one declining faster. When the pressure becomes sufficiently anisotropic, with P_∥,e > P_⊥,e , the plasma becomes unstable to the SFHI (we confirmed separately that the instability is absent for simulations of comparable durations if the initial β is too small, β₀ ≲ 1). The evolution of the volume-averaged electron pressure ratio 〈P_⊥,e /P_∥,e 〉 versus the volume-averaged electron plasma beta 〈β_e,∥〉 is shown in the top panel of Figure 3 for all four compositions. In the three charged cases, the evolution of 〈P_⊥,e /P_∥,e 〉 traces the expected early-time scaling P_⊥,e /P_∥,e ∼ 1 − (8/15)(t/τ_cool) (Equation (8)) until the onset of the SFHI. Thereafter, the decrease of 〈P_⊥,e /P_∥,e 〉 stalls while β_e,∥ continues to decline, and the trajectory fluctuates along the corresponding instability threshold. The pair-plasma case lies near the instability threshold of Equation (9) with C_thr ≈ 0.7 and the semi-relativistic electron–ion case near C_thr ≈ 1.4, in agreement with the analytical expectations from Section 2.2. Meanwhile, the relativistic electron–ion case is closest to the threshold with C_thr ≈ 2, somewhat larger than the expectation of C_thr ≈ 1.4; one reason for this discrepancy may be the influence of thermal coupling between the electrons and relativistic ions, which is discussed in Section 3.3.

Zoom In Zoom Out Reset image size

In general, the firehose instability can be either parallel or oblique, with the dominant orientation determined by the physical parameters and plasma composition. In our PIC simulations, we find that the properties of the magnetic-field fluctuations produced by the SFHI depend on the plasma composition. The profiles of the perpendicular components B_x and B_z are shown in Figure 4 for the three unstable compositions at double the fiducial size, taken at representative times after the onset of the instability (fluctuations in B_y are negligible). The pair-plasma case (top row) has a mixture of parallel and oblique modes; the parallel mode is visible in the B_x profile while the oblique mode distorts the B_z profile. The relativistic electron–ion case (middle row) is dominated by circularly polarized parallel fluctuations, with no visible signatures of oblique modes. The case with sub-relativistic ions (bottom row), on the other hand, is strongly dominated by oblique modes (in B_z ), with no discernible sign of parallel modes (in B_x ); this is consistent with expectations from the nonrelativistic regime, in which the oblique firehose generally grows faster than the parallel firehose (e.g., Li & Habbal 2000; Bott et al. 2021).

Zoom In Zoom Out Reset image size

The properties of the fluctuations also change over time: they typically start at a wavelength that is approximately twice larger than 2π ρ_e0, then distort and cascade to smaller-scale waves. In Figure 5, we show the evolution of the magnetic energy spectrum (in the k_y direction, integrated over k_x and averaged over five time snapshots) at early times for the medium-size pair-plasma case (L/2π ρ_e0 ≈ 8.1). The initial state (a flat spectrum) is due to numerical PIC noise. The early-time peak, during the linear stage, is at a wavenumber of k_y ≈ 0.4/ρ_e0, broadly consistent with previous investigations of the kinetic electron-firehose instability (e.g., Riquelme et al. 2018; Kim et al. 2020; López et al. 2022). Note that the linear theory from the Appendix cannot be used to predict the most unstable wavenumber, due to the long-wavelength/low-frequency approximation. The peak then moves to higher wavenumbers and decays (the longer-term evolution of the system is described further below). There does not appear to be a significant inverse transfer of energy to larger scales. The evolution of the spectrum is qualitatively similar for all of the electron–ion cases, peaking at a similar value of k ρ_e0, which confirms the instability as being the electron firehose. In principle, since the electron pressure anisotropy causes an anisotropy in the total plasma pressure, it may trigger the fluid firehose instability, which would occur at ion scales of k_y ≲ 1/ρ_i rather than the electron scales. However, this is not observed in our PIC simulations (even in the large-domain case having L/2π ρ_i0 ≈ 3.5), presumably because the electron-firehose instability onset and growth is more rapid, and therefore depletes the pressure anisotropy before the ion-scale fluid firehose modes can grow. This picture is supported by recent work that considers coexisting electron and ion firehose instabilities in nonrelativistic plasmas, which finds that electron-firehose modes grow faster than ion-scale ones (López et al. 2022).

Zoom In Zoom Out Reset image size

The onset of the SFHI is associated with a rapid increase in the magnetic energy of the perturbations (from the initial noise floor), as shown in the top panel of Figure 6 for all three charged cases. The instability starts to grow roughly at the time τ_onset predicted by the linear theory, except for the semi-relativistic case, where there is a substantial delay (possibly indicating an overshoot, thus requiring smaller τ_gyro/τ_cool to get into the asymptotic regime).

Zoom In Zoom Out Reset image size

The initial SFHI growth can be studied in greater detail by subtracting off the contribution to magnetic energy from the PIC noise (caused by the finite number of particles per cell), which we denote E_noise. In each case, we define E_noise as the mean magnetic energy during an interval of duration 0.015τ_cool shortly before the instability starts to grow. Physically, E_noise is associated with fluctuations over a broad range of wavenumbers, as seen from the flat spectrum during the earliest (pre-instability) time in Figure 5; E_noise overwhelms the early-time signal of the instability in E_mag, even if the instability dominates the energy at the most unstable wavenumber. Subtracting E_noise from E_mag allows the initial SFHI growth to be measured over several orders of magnitude in energy, as shown in the middle panel of Figure 6 (note that this plot shows the medium-size simulations for clarity). We find that this growth can be represented very well by either an exponential or a super-exponential function, as predicted by the linear theory in Section 2.3. Motivated by Equation (13) describing unregulated growth, we choose to fit the data with a function of the form

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\mathrm{mag}}(t)}{{E}_{\mathrm{mag},\mathrm{mean}}}\sim \eta \exp \left[{\left(\displaystyle \frac{t-{\tau }_{0}}{{\tau }_{1}}\right)}^{3/2}\right],\end{eqnarray} \tag{ 14 }$$

where η, τ₀, and τ₁ are fitting parameters representing the amplitude, onset time, and growth time, respectively. Fits of this form are shown for each case by black dashed/dotted lines in the middle panel of Figure 6. For all cases, we choose η = 2 × 10⁻⁶, representing the level of the seed fluctuations from the noise. Fastest onset and growth occurs for the pair-plasma case, with τ₀ = 0.065τ_cool ≈ τ_onset,pair (where τ_onset,pair is the theoretical linear onset time from Equation (11) with C_thr = 0.7 and β_e0 = 20) and τ₁ ≈ 9 × 10⁻³ τ_cool. This is in close agreement with the linear theory (Equations (12) and (13)), which predicts that τ₁ = 2^−2/3 τ_gr = 7 × 10⁻³ τ_cool for the given plasma parameters and assuming k ρ_e0 = 0.4. The next fastest growth is for the relativistic electron–ion case, with τ₀ = 0.11τ_cool ≈ 0.85τ_onset,ei (where τ_onset,ei is the linear onset time with C_thr = 1.4 and β_e0 = 20) and τ₁ = 1.6 × 10⁻² τ_cool, such that the growth rate is nearly a factor of 2 lower than in the pair case. For comparison, the linear theory calculation predicts τ₁ = 9 × 10⁻³ τ_cool in this case (again assuming k ρ_e0 = 0.4), somewhat faster than measured. Finally, the semi-relativistic plasma has the slowest growth, with τ₀ = 0.21τ_cool ≈ 1.6τ_onset,ei and τ₁ ≈ 2.6 × 10⁻² τ_cool; incidentally, this growth rate is approximately a factor of 3 lower than the pair case. The growth rate is significantly below the linear prediction, τ₁ = 1.8 × 10⁻² τ_cool. In summary, the hierarchy of the growth rates is qualitatively consistent with the estimates from linear theory, with pair plasma having the closest agreement with theory. Deviations from linear theory for the electron–ion cases may be due to the non-asymptotic nature of the simulations, as well as corrections from higher-order kinetic terms and the oblique nature of the dominant modes.

Following the initial super-exponential growth, the magnetic energy saturates and then decays. On a much longer timescale, however, the magnetic energy undergoes a nonlinear cyclical evolution of growth and decay (as seen in the top panel of Figure 6). The period of the cycles is fastest for the pair-plasma case (having a duration of τ_cycle ≈ τ_cool/12) and slowest for the case with sub-relativistic ions (with τ_cycle ≈ τ_cool/3). Each burst has a peak magnetic energy that is a small fraction of the energy in the mean magnetic field, ≲0.02E_mag,mean, and rises above the local minima (between peaks) by a factor of ∼3. The peak amplitudes generally decay over an even longer timescale (of order ∼τ_cool). The peak amplitude appears to be sensitive to parameters that are not the focus of our study, such as τ_cool/τ_gyro (not shown). Eventually, because of the continuous decrease of β_∥,e from cooling, the plasma will exit the firehose unstable region of parameter space (and presumably return to a stable cooling evolution).

In principle, the SFHI may affect the cooling rate of the plasma. Pitch-angle scattering will mix the weakly cooled, parallel-propagating particles with strongly cooled, perpendicular-propagating particles, thereby enhancing the cooling rate over the predicted stable evolution (additionally, some electron kinetic energy is transferred to firehose magnetic fields, but this has a weaker impact on the cooling rate). An enhancement of the cooling rate (over the uncharged evolution) is indeed observed for the pair-plasma and semi-relativistic electron–ion cases, as shown in the bottom panel of Figure 6, which displays the evolution of the total kinetic energy E_kin,s for each species (s ∈ {e, i}). The relativistic electron–ion case, on the other hand, experiences slower cooling of the electrons than the uncharged simulation. This implies that electrons experience anomalous heating from some source. The only available source that electrons can draw energy from is the ion kinetic energy; indeed, we find that the anomalous electron heating is compensated by a decrease of ion kinetic energy (dashed blue line in Figure 6), indicating collisionless thermal coupling between the two species (studied in Section 3.3). The adjustment to the cooling rate caused by the SFHI in all cases is small (approximately ∼10% over the simulated timescales), but in principle may build to a significant temperature difference over time. We also note that there does not appear to be a simple analytical function that describes E_kin,s (t): power laws and exponential functions provide poor fits to the simulation data.

The electron momentum distribution arising from the SFHI is close to an isotropic thermal distribution, as shown separately for the parallel and perpendicular components in Figure 7 (red lines) for the pair-plasma case at t/τ_cool = 1.4 (other cases are similar). The stable uncharged simulation, by comparison, has a very strongly nonthermal distribution at this time (blue lines in Figure 7, which closely follow the analytical solution from Figure 1). Thus, the SFHI efficiently thermalizes the distribution by mixing particles with different pitch angles. At large momenta, there is a modest excess of particles in the parallel direction compared to the perpendicular direction, highlighting a persistent asymmetry caused by the continuous synchrotron cooling, of the amount necessary to maintain pressure anisotropy at the marginally unstable state. Although the overall distribution is fairly close to an isotropic Maxwell–Jüttner distribution, the anisotropy is not represented well by a bi-Maxwell–Jüttner distribution (i.e., the anisotropies are comparable to the deviations from the Maxwell–Jüttner distribution in each direction, indicating that the anisotropic component is nonthermal).

Zoom In Zoom Out Reset image size

We next study the effect of the SFHI on ions. In the pair-plasma case, the positrons undergo an evolution that is identical to that of the electrons (as required by symmetry), which was already described in Section 3.2. Therefore, in this subsection, we consider what happens in the simulations with non-radiative ions.

The ions remain in their initial thermal state until the SFHI produces magnetic fields that can perturb them. For the ultra-relativistic electron–ion case, some time after the onset of SFHI, the ions grow a positive pressure anisotropy, P_⊥,i > P_∥,i , as shown in the bottom panel of Figure 3. This occurs because the ions interact with the fields produced by the electron firehose.

The initial growth of the ion pressure anisotropy can be modeled in the limit of double-adiabatic evolution in an increasing magnetic field. In this simplified model, the ions conserve their magnetic moment $\mu ={p}_{\perp }^{2}/2{m}_{i}B$ (associated with periodic Larmor orbits) and parallel momentum p_∥ (associated with periodic motion along the magnetic-field lines) as the firehose fluctuations increase the magnetic-field strength. Under this "double-adiabatic" evolution, the ion distribution will evolve as

$$\begin{eqnarray}&&{f}_{i,\mathrm{adi}}({p}_{\perp },{p}_{\parallel },t)=\displaystyle \frac{{B}_{0}}{B(t)}{f}_{i0}\left(\sqrt{\displaystyle \frac{{B}_{0}}{B(t)}{p}_{\perp }^{2}+{p}_{\parallel }^{2}}\right),\end{eqnarray} \tag{ 15 }$$

where B(t) is the orbit-averaged magnetic field (which we approximate as being the same for all particles) and f_i0(p) is the initial isotropic distribution. Assuming a weak fluctuation, B/B₀ = 1 + _{δ B} with _{δ B} ≪ 1, we can expand the distribution in Equation (15) and compute the pressures to show that the adiabatic evolution follows

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{P}_{\perp ,i}}{{P}_{\parallel ,i}}\right|}_{\mathrm{adi}}=1+\displaystyle \frac{2}{5}{\epsilon }_{\delta B}=3-2\displaystyle \frac{{\beta }_{\parallel ,i}}{{\beta }_{\parallel ,i,0}}.\end{eqnarray} \tag{ 16 }$$

In deriving Equation (16), we have assumed that the ions are ultra-relativistic. The double-adiabatic evolution given by Equation (16) is shown by the cyan line in the bottom panel of Figure 3. The initial evolution of the simulation is close to this adiabatic prediction, but shifted slightly.

Departures from the double-adiabatic scaling are due to pitch-angle scattering events that break the conservation of μ and p_∥. To demonstrate this, in Figure 8 we show the evolution of the magnetic moment μ(t) for 11 randomly selected ions (from the ultra-relativistic electron–ion case). To reduce the effect of noise, these ions are chosen with the criterion that their initial momentum is larger than twice the average ion momentum, $p(t=0)\gt 2{\overline{p}}_{0}$. We find that μ is well conserved up to a time of t/τ_cool ∼ 0.15 (with small fluctuations due to noise), at which point electron-firehose fluctuations emerge. Subsequently, the conservation of μ is broken for most particles. Note that high-frequency oscillations in μ are observed for some particles due to their Larmor gyrations; the true adiabatic invariant should be evaluated at the guiding-center position rather than particle position, and would presumably exhibit a smoother evolution.

Zoom In Zoom Out Reset image size

Returning to the bottom panel of Figure 3, once the ion anisotropy becomes sufficiently large, either the ion mirror instability or ion-cyclotron instability may be triggered. The threshold for the (nonrelativistic) mirror instability with simultaneous anisotropy in both electrons and ions can be written as (Stix 1962; Hellinger 2007)

$$\begin{eqnarray}&&\begin{array}{l}{\beta }_{\parallel ,e}\displaystyle \frac{{P}_{\perp ,e}}{{P}_{\parallel ,e}}\left(\displaystyle \frac{{P}_{\perp ,e}}{{P}_{\parallel ,e}}-1\right)+{\beta }_{\parallel ,i}\displaystyle \frac{{P}_{\perp ,i}}{{P}_{\parallel ,i}}\left(\displaystyle \frac{{P}_{\perp ,i}}{{P}_{\parallel ,i}}-1\right)\\ \quad \gt 1+\displaystyle \frac{{\left({P}_{\perp ,i}/{P}_{\parallel ,i}-{P}_{\perp ,e}/{P}_{\parallel ,e}\right)}^{2}}{2\left(1/{\beta }_{\parallel ,e}+1/{\beta }_{\parallel ,i}\right)},\end{array}\end{eqnarray} \tag{ 17 }$$

while the one for the ion-cyclotron instability is typically fit by ${P}_{\perp ,i}/{P}_{\parallel ,i}\gtrsim 1+A/{\beta }_{\parallel ,i}^{b}$, with A a free coefficient and b ≈ 0.4 (e.g., Gary & Lee 1994; Hellinger et al. 2006; Sironi & Narayan 2015). These thresholds do not account for relativistic effects, departures from bi-Maxwellian distributions, and reduced effective (relativistic) mass ratios, but we note that López et al. (2016) found that similar thresholds apply for the cyclotron instability in relativistic pair plasmas with anisotropy in a single species, which resembles our setup here. We show these thresholds in the bottom panel of Figure 3. For the mirror threshold, we choose representative electron parameters of P_⊥,e /P_∥,e = 0.9 and β_∥,e = 0.18 in Equation (17); this threshold is not reached by the simulation. For the ion-cyclotron threshold, we choose A = 0.31; the simulation crosses this threshold and then alters its trajectory (times of t/τ_cool = 0.4 and t/τ_cool = 0.8 are marked, with the latter corresponding to the time at which the ion pressure anisotropy becomes limited). This suggests that the ion-cyclotron instability is the dominant secondary instability in our simulations. Further evidence for the lack of efficient mirror growth comes from the polarization of the modes: we do not observe parallel magnetic-field fluctuations (δ B_y ) or density fluctuations that would be expected from the mirror instability. Rather, the polarization of the modes remains similar to the electron-firehose modes, as expected for the ion-cyclotron instability.

The above suggests that the electron–ion thermal coupling observed in our PIC simulations is caused by the ion-cyclotron instability. Free energy for the ion-cyclotron magnetic-field fluctuations is provided by the ions' internal energy (in the form of ion pressure anisotropy). The dissipation of these fluctuations will cause the heating of electrons and ions (partitioned in an undetermined way). The end result of this process is the net heating of electrons and cooling of ions.

Throughout the evolution, the ion momentum distribution remains close to a thermal distribution in each direction. This is in contrast with the moderate nonthermal anisotropy for the radiatively cooled electrons.

While ions also attain P_⊥,i /P_∥,i > 1 in the semi-relativistic case, the pressure anisotropy does not cross the ion instability thresholds and thus the amount of thermal coupling is negligible in this case. The limited amount of ion pressure anisotropy in this case may be a consequence of the scale mismatch between the ion gyroradius ρ_i and electron gyroradius ρ_e , having ${\rho }_{e}/{\rho }_{i}\sim ({T}_{e}/{T}_{i}){\theta }_{i}^{1/2}$ in the semi-relativistic regime. When ρ_i ≫ ρ_e , the ions are unable to interact strongly with the electron-scale firehose modes (in contrast to the relativistic ion case where ρ_i ∼ ρ_e ). Thus, we conclude that electron–ion thermal coupling is limited to the relativistic ion regime in which ρ_i ∼ ρ_e .