A nonlinear Monte Carlo model of super-diffusive shock acceleration with magnetic field amplification

TL;DR

This paper introduces a nonlinear Monte Carlo model for super-diffusive shock acceleration that captures magnetic field amplification and anisotropic particle distributions, providing insights into shock precursor dynamics and magnetic turbulence.

Contribution

It presents a novel Monte Carlo approach that models super-diffusive particle transport and magnetic field amplification in collisionless shocks, including full angular anisotropy.

Findings

Super-diffusive transport causes strong quadruple anisotropy in particle distribution.

Pressure anisotropy leads to non-resonant mirror-type instability.

Model offers a comprehensive description of shock structure and magnetic turbulence.

Abstract

Fast collisionless shocks in cosmic plasmas convert their kinetic energy flow into the hot downstream thermal plasma with a substantial fraction of energy going into a broad spectrum of superthermal charged particles and magnetic fluctuations. The superthermal particles can penetrate into the shock upstream region producing an extended shock precursor. The cold upstream plasma flow is decelerated by the force provided by the superthermal particle pressure gradient. In high Mach number collisionless shocks, efficient particle acceleration is likely coupled with turbulent magnetic field amplification (MFA) generated by the anisotropic distribution of accelerated particles. This anisotropy is determined by the fast particle transport making the problem strongly nonlinear and multi-scale. Here, we present a nonlinear Monte Carlo model of collisionless shock structure with super-diffusive propagation of high-energy Fermi accelerated particles coupled to particle acceleration and MFA which affords a consistent description of strong shocks. A distinctive feature of the Monte Carlo technique is that it includes the full angular anisotropy of the particle distribution at all precursor positions. The model reveals that the super-diffusive transport of energetic particles (i.e., Levy-walk propagation) generates a strong quadruple anisotropy in the precursor particle distribution. The resultant pressure anisotropy of the high-energy particles produces a non-resonant mirror-type instability which amplifies compressible wave modes with wavelengths longer than the gyroradii of the highest energy protons produced by the shock.

Equations

$$⟨\Delta z^2⟩=A_D{t}^b,$$

$$\rho (z)u(z)={\rho }_0{u}_0,$$

$${\Phi }_P^{\mathrm{part}}(z)+P_w(z)={\Phi }_{P0},$$

$$\rho (z)u^2(z)+P_{\mathrm{th}}(z)+P_{\mathrm{ep}}(z)+P_w(z)={\Phi }_{P0},$$

$${\Phi }_E^{\mathrm{part}}(z)+F_w(z)={\Phi }_{E0},$$

$$\frac{\rho (z)u^3(z)}{2}+F_{\mathrm{th}}(z)+F_{\mathrm{ep}}(z)+F_w(z)+Q_{\mathrm{esc}}={\Phi }_{E0},$$

$$F_{\mathrm{th}}(z)=u(z)\frac{{\gamma }_g{P}_{\mathrm{th}}(z)}{{\gamma }_g-1},$$

$$u(z)\frac{\partial W\left(z,k\right)}{\partial z}+\frac{3}{2}W\left(z,k\right)\frac{du(z)}{dz}+\frac{\partial \Pi \left(z,k\right)}{\partial k}=$$

$$=G(z,k)-\mathcal{L}(z,k),$$

$$\Pi (z,k)=-\frac{D_{\mathrm{K}}}{\sqrt{\rho (z)}}{k}^{\frac{11}{2}}W(z,k{)}^{\frac{1}{2}}\frac{\partial }{\partial k}\left[\frac{W(z,k)}{k^2}\right],$$

$${\int }_{k_{\min }}^{k_{\max }}W(z_{\mathrm{FEB}},k)dk=\frac{B_0^2}{4\pi }.$$

$$\delta {\theta }_{\max }=\sqrt{6\delta t/t_c},$$

$$\Psi ({\lambda }_{\mathrm{LF}})\propto {\lambda }_{\mathrm{LF}}^{-\nu }\text{for}{\lambda }_{\mathrm{LF}}(z,p)>{\lambda }_0(z,p).$$

$${\lambda }_{\mathrm{LF}}(z,p)=\tan \left(\frac{\pi {\xi }_0}{2}\right){\lambda }_0(z,p)\mathrm{for}z<z_{\mathrm{LF}}.$$

$$F({\lambda }_{\mathrm{LF}})=\frac{2}{\pi }\arctan \left(\frac{{\lambda }_{\mathrm{LF}}}{{\lambda }_0}\right),$$

$$\Psi ({\lambda }_{\mathrm{LF}})=\frac{2}{\pi {\lambda }_0}{\left[1+{\left(\frac{{\lambda }_{\mathrm{LF}}}{{\lambda }_0}\right)}^2\right]}^{-1}.$$

$${\lambda }_{\mathrm{LF}}={\lambda }_0\left[\frac{2\nu \left(\nu -2\right){\xi }_1}{{\left(\nu -1\right)}^2}+1\right],{\xi }_1\le D_0,$$

$${\lambda }_{\mathrm{LF}}={\lambda }_0\left[\frac{2(\nu -2)}{\nu -1}{\left[\frac{1}{\nu \left(1-{\xi }_1\right)}\right]}^{\frac{1}{\nu -1}}+1\right],{\xi }_1>D_0,$$

$$\Psi \left({\lambda }_{\mathrm{LF}}^{\ast }\right)=\frac{C_{\nu }}{{\lambda }_0},\mathrm{for}{\lambda }_{\mathrm{LF}}^{\ast }\le {\lambda }_{\ast },$$

$$\Psi \left({\lambda }_{\mathrm{LF}}\right)=\frac{C_{\nu }}{{\lambda }_0}{\left(\frac{{\lambda }_{\mathrm{LF}}^{\ast }}{{\lambda }_{\ast }}\right)}^{-\nu },\mathrm{for}{\lambda }_{\mathrm{LF}}^{\ast }>{\lambda }_{\ast },$$

$${\lambda }_{\ast }=2{\lambda }_0\frac{\nu -2}{\nu -1}\text{and}{C}_{\nu }=\frac{(\nu -1)D_0}{2(\nu -2)}$$

$${\int }_0^{\infty }\Psi \left(\lambda \right)d\lambda =1$$

$${\int }_0^{\infty }\lambda \Psi \left(\lambda \right)d\lambda ={\lambda }_0,$$

$$\delta t(z,p)=\frac{A_{\theta }^2{\lambda }_{\mathrm{LF}}(z,p)}{6v_{\mathrm{pf}}},$$

$$Q_{\mathrm{esc}}(p)=-\frac{J(z_{\mathrm{FEB}},p)p^4}{4\pi m_{p}c{u}_0},$$

$$J(z,p)=2\pi {\int }_{-1}^{1}d\mu v\mu f_{\mathrm{pf}}^{\mathrm{ep}}(z,p,\mu ),$$

$$A_1(z)=N^{-1}{\int }_0^{\infty }J(z,p)v^{-1}{p}^{2}dp.$$

$$A_2(z,p)=\pi {\int }_{-1}^1(3{\mu }^2-1)f_{\mathrm{pf}}^{\mathrm{ep}}(z,p,\mu )d\mu ,$$

$$A_2(z)=N^{-1}{\int }_0^{\infty }{A}_2(z,p)p^{2}dp,$$

$$N=2\pi {\int }_0^{\infty }{p}^{2}dp{\int }_{-1}^{1}d\mu f_{\mathrm{pf}}^{\mathrm{ep}}(z,p,\mu ).$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\nabla )\mathbf{u}\right)=-\nabla p_g+\frac{1}{4\pi }(\nabla \times \mathbf{B})\times \mathbf{B}-$$

$$-\frac{1}{c}({\mathbf{j}}^{\mathrm{ep}}-e{n}_{\mathrm{ep}}\mathbf{u})\times \mathbf{B}-\int \mathbf{p}I[f^{\mathrm{ep}}]d^{3}p,$$

Full text

\pacs

PACS: \pacsValid PACS appear here

1. Introduction

In contrast to collision-dominated shocks, strong collisionless plasma shocks are capable of converting the kinetic power of the upstream flow to both thermal and non-thermal components. The conversion process, in diffuse plasmas where Coulomb collisions are very infrequent, is due to highly nonlinear interactions between the particles and the background magnetic turbulence utilizing the first-order Fermi mechanism (Balogh & Treumann, 2013; Marcowith et al., 2016). While collisionless shocks are difficult to study in the laboratory, they are known to exist in cosmic plasmas and play a critical role in producing nonthermal particles observed throughout the cosmos. The lack of Coulomb collisions allows the nonthermal components – accelerated energetic particles (EP) and magnetic turbulence – to be long-lived and dynamically significant. Observations and theory both confirm that particle acceleration can be efficient enough so the energetic particles that penetrate into the shock precursor can slow the bulk supersonic flow significantly before the viscous subshock occurs. The viscous subshock is a small-length-scale collisionless shock of moderate Mach number $\sim 3$. The subshock, which is directly observed in heliospheric shocks and particle-in-cell (PIC) simulations, involves mainly thermal particles and is required to produce the entropy and jumps in plasma density, temperature, and magnetic field needed to satisfy the Rankine-Hugoniot relations. The superthermal precursor particles are highly anisotropic and drive a number of plasma instabilities (e.g., Bell, 2004; Schure et al., 2012; Bykov et al., 2013) producing magnetic turbulence which interacts with the precursor particles producing the turbulence. Modeling strong collisionless shocks is an interesting multi-scale problem strictly from the point of view of nonlinear plasma physics. It is also fundamentally important for understanding energetic particle populations observed in space.

The fact that collisionless shocks accelerate ambient particles in many locations ranging from the Earth bow shock to shocks in galactic clusters is widely accepted. While a great deal is known about the acceleration process, the most important part—the collisionless wave-particle interactions driving particle isotropization—remains uncertain. In principle, a full description of the plasma interactions is obtainable with PIC simulations (e.g., Park et al., 2015). However, PIC simulations are computationally expensive and results thus far are limited to a relatively narrow dynamical range which is particularly restricting for the modeling of nonrelativistic  shocks such as those seen in supernova remnants.

All collisionless shock calculations, other than PIC, must approximate particle transport and most models assume fast particles obey standard diffusion (in the local plasma rest frame) where the mean-square-displacement is proportional to time, i.e.,

[TABLE]

with $b=1$. This simple equation (even with $b=1$) hides a great deal of complexity since the proportionality factor, $A_{D}$, depends non-linearly on the details of the self-generated magnetic turbulence and will vary with particle momentum, position relative to the subshock, and the Fermi acceleration efficiency.

Furthermore, there is no fundamental reason why $b=1$ in Eq. (1) and both super- ($b>1$) and sub-diffusive ($b<1$) transport regimes are possible in complex nonlinear and intermittent systems (e.g., Shlesinger et al., 1993; Zumofen & Klafter, 1993; Zelenyi & Milovanov, 2004; Zaburdaev et al., 2015). Such non-standard diffusion has been shown to be important in laboratory and fusion plasmas (e.g., Balescu, 2005; Perrone et al., 2013; Bovet et al., 2015). Furthermore, there is evidence from in-situ spacecraft observations of heliospheric shocks for anomalous diffusion where the mean-square-displacement grows non-linearly with time with $b>1$ (Perri et al., 2015). The effect of complex transport on EP propagation and acceleration was discussed by (Kirk et al., 1996; Perri & Zimbardo, 2009; Malkov & Diamond, 2009; Zimbardo & Perri, 2013; Lazarian & Yan, 2014).

We consider Fermi acceleration in strong quasi-parallel shocks where the average magnetic field direction upstream of the shock is close to the shock normal, implicitly assuming the magnetic field at the subshock is turbulent enough so effects from perpendicular components of the field can be ignored (e.g., Ellison et al., 1995). An analysis of obliquity effects on particle propagation within the Monte Carlo model suggests that shocks can be considered “parallel” for angles up to $\sim\pi/4$ from the shock normal. The shock produces anisotropic EP distributions in the shock precursor that result in strong non-adiabatic amplification of irregular magnetic fields by EP-driven instabilities (e.g., Marcowith et al., 2016). Magnetic field fluctuations present in the interstellar medium are highly amplified by these instabilities as they traverse the shock precursor.

An important characteristic of strong shocks undergoing efficient Fermi acceleration is that the highest energy particles are distributed throughout the entire precursor while lower energy particles are concentrated close to the subshock. This results in a strong spatial dependence of the growth rate of magnetic fluctuations with a given wavenumber and may lead to super-diffusion in the outer precursor where the magnetic turbulence is growing and is likely highly intermittent.

Here, we assume the intermittency of the turbulence dominates the EP propagation in the super-diffusion region of the precursor before the growing turbulence saturates. The boundary of the super-diffusion region in the upstream is parameterized by $z_{\mathrm{LF}}$ which is the distance from the subshock. This parameterization is needed since we use simplified models for magnetic turbulence cascade, which deal only with the spectrum of turbulence, while the intermittency requires a more detailed description. This is discussed in more detail in Section 2.

Super-diffusive EP propagation has been seen at a few gyro-rotation periods in magnetohydrodynamic (MHD) simulations of the non-resonant EP-driven instability in (Reville et al., 2008) where the EP current was fixed (i.e., without considering self-consistent EP evolution). However, the question whether EP transport is diffusive in the shock precursor when EP-driven instabilities are rapidly growing has not yet been addressed. The nonlinear backreaction of EPs on the shock structure will influence the EP transport and anisotropy making nonlinear calculations coupling EP production, shock structure, instability growth, and particle transport essential. While in principle, the full time-dependent picture can be modeled with large-scale PIC simulations, such calculations are still beyond current capabilities. Until these calculations become available, the effect of non-diffusive propagation on nonlinear Fermi shock acceleration can be studied with Monte Carlo simulations of particle transport by simply assuming that the EP transport in some regions of the shock precursor is non-diffusive.

An important distinction with Fermi shock acceleration models which are based on the advection-diffusion transport equation is that they model only EP-current anisotropies. The Monte Carlo model does not make a diffusion approximation and thus can account for arbitrary angular anisotropy harmonics. We find that the inclusion of super-diffusion from Lévy-walk scattering-length distributions produces specific anisotropies in the particle distributions that drive instabilities that do not occur with standard diffusion models.

When Fermi shock acceleration is efficient, particle transport, the shock structure, MFA, turbulence cascading, and thermal particle injection must all be calculated self-consistently. No technique, not even PIC simulations, can currently do this full calculation from first principles over a dynamic range sufficient to model EP production in SNRs or other strong nonrelativistic shocks. While approximations must be made, much of the essential nonlinear physics can be modeled with Monte Carlo techniques.

The Monte Carlo simulation we use couples analytic descriptions of resonant and non-resonant wave growth with anisotropic particle transport in EP-dominated shocks (see Ellison et al., 1996; Bykov et al., 2014, and references therein). All of the nonlinear effects mentioned above have been consistently included assuming standard diffusion. We now generalize the Monte Carlo model by explicitly including super-diffusion in the shock precursor. Because the Monte Carlo model accounts for the full anisotropic EP distribution functions, the dispersion relations we derive simultaneously include the EP-driven resonant streaming instability, and the two EP-current driven instabilities: Bell’s short-wavelength instability, and the long-wavelength instability (see (Schure et al., 2012; Bykov et al., 2013)).

Our nonlinear model shows two distinctive features. The first is that super-diffusion results in a highly amplified specific quadrupole anisotropy of EP particles. This anisotropy produces a mirror instability that has not been previously considered in efficient shock acceleration. The mirror instability contributes significantly to the generation of long-wavelength magnetic turbulence which, in principle, can be studied with direct measurements of heliospheric shocks and with indirect analysis of the broadband synchrotron emission seen in supernova remnants. The second feature is that super-diffusion in the shock precursor results in a substantial broadening of the spectrum of energetic particles escaping the precursor. Energetic particles escaping the shock will undergo inelastic collisions with surrounding matter and produce high-energy radiation. We quantify the mirror instability and escaping EPs with a limited number of Monte Carlo examples.

2. The nonlinear Monte Carlo shock model

We construct a steady-state model of a plane-parallel, nonrelativistic collisionless shock where the nonlinear shock structure is determined iteratively. The shock is directed in the $-z$-direction with a subshock at $z=0$ and an upstream free escape boundary (FEB) limiting Fermi acceleration at $z=z_{\mathrm{FEB}}$. For simplicity, we model only protons since they mainly determine the nonlinear shock structure and drive the long-wavelength instabilities we consider. Electrons can be included, as in Warren et al. (2015), when radiation is calculated.

The Monte Carlo shock model includes the following main elements:

(1) Particle injection, which is self-consistently coupled to Fermi acceleration where some fraction of shock-heated thermal particles re-cross the subshock (assumed transparent), gain additional energy, and enter the acceleration process. Any particle that crosses from downstream back upstream at least once is termed energetic;

(2) Shock-smoothing, where backpressure from superthermal particles slows and heats the precursor plasma upstream of the viscous subshock in order to conserve momentum and energy;

(3) The self-consistent  determination of the overall shock compression ratio, $R_{\mathrm{tot}}$, taking into account escaping EPs, magnetic pressure, and the modification of the equation of state from the production of relativistic particles;

(4) Fluctuating magnetic fields simultaneously calculated from resonant, short-wavelength, long-wavelength, and mirror instabilities generated from the EP current, and super-diffusion pressure anisotropies in the shock precursor;

(5) Momentum and position dependent particle transport determined from the self-generated magnetic turbulence;

(6) A determination of the local plasma scattering center speed relative to the bulk plasma from energy conservation without assuming Alfvén waves; and,

(7) Turbulence convection and compression, cascade, and dissipation of wave energy into the background plasma. The iterative Monte Carlo approach allows all of these processes to be coupled and calculated simultaneously in a reasonably consistent fashion.

3. Mass-Energy-Momentum Conservation

We determine the self-consistent shock structure with an iterative procedure by forcing mass-energy-momentum conservation. All particles–thermal and superthermal–are transported through the shock keeping full account of the anisotropic particle distribution and the momentum and energy contributions from the magnetic fluctuations (see Bykov et al. (2014) for full details).

In the shock rest frame, the mass flow conservation is given by

[TABLE]

where $\rho(z)$ is the plasma density, $u(z)$ is the bulk flow speed, and the subscript “0” here and elsewhere indicates far upstream values. We define the “shock structure” as $u(z)$, where $z$ is the distance measured from the subshock at $z=0$. The momentum flux conservation is determined by

[TABLE]

where $\Phi^{\mathrm{part}}_{P}(z)$ is the particle momentum flux, $P_{w}(z)$ is the momentum flux carried by the magnetic turbulence, and $\Phi_{P0}$ is the far upstream momentum flux, i.e., upstream from the free escape boundary where the interstellar magnetic field is $B_{0}$.

Separating the contributions from the thermal and accelerated particles we have

[TABLE]

where $P_{\mathrm{th}}(z)$ is the thermal particle pressure and $P_{\mathrm{ep}}(z)$ is the accelerated particle pressure. As mentioned above, a particle is “accelerated” if it has crossed the subshock more then once. There is no other injection threshold and even though we use the subscript “EP”, the vast majority of accelerated particles will always be nonrelativistic. Of course, if the acceleration is efficient, a large fraction of the pressure may be in relativistic particles.

The energy flux conservation law is

[TABLE]

where $\Phi^{\mathrm{part}}_{E}(z)$ and $F_{w}(z)$ are the energy fluxes in particles and magnetic field correspondingly, and $\Phi_{E0}$ is the energy flux far upstream. Taking into account particle escape at an upstream FEB, this can be re-written as

[TABLE]

where $F_{\mathrm{th}}(z)$ is the internal energy flux of the background plasma, $F_{\mathrm{ep}}(z)$ is the energy flux of accelerated particles, and $Q_{\mathrm{esc}}$ is the energy flux of particles that escape at the upstream FEB (note that $Q_{\mathrm{esc}}$ is defined as positive even though EPs escape moving in the negative $z$-direction).

The separation between “thermal” particles and “accelerated” particles in a shock undergoing diffusive shock acceleration is not necessarily well defined. Furthermore, the energy exchange between the thermal and superthermal populations is certain to occur through non-trivial wave-particle interactions. Nevertheless, the bulk of the plasma mass will always be in quasi-thermal background particles and the internal energy flux of this background plasma can be expressed as

[TABLE]

where $\gamma_{g}=5/3$ is the adiabatic index of the background thermal plasma.

All of the quantities in Eqs. (2)–(7) are directly measured in the Monte Carlo simulation. The magnetic turbulence, and therefore $P_{w}(z)$ and $F_{w}(z)$, is determined from $u(z)$ and the analytic expressions for wave growth and cascading discussed below. Once the assumptions for wave growth are made, the equations for mass, momentum, and energy flux are over determined and a unique, nonlinear solution conserving mass, momentum, and energy flux in the shock rest frame can be found by iterating $u(z)$. In practice, a “consistent solution” is accepted when the momentum and energy fluxes are conserved to within a few percent at all $z$.

3.1. Turbulence cascade

The magnetic turbulence energy flow $F_{w}(z)$ in Eq. (6) is determined by the spectral energy density of the magnetic fluctuations $W(z,k)$ (see Bykov et al. (2014) for details ) which obeys the equation

[TABLE]

where $\Pi(z,k)$ is the flux of magnetic energy, $G(z,k)$ is the spectral energy growth rate due to EP instabilities, and $\mathcal{L}(z,k)$ is the turbulence dissipation rate. Following Matthaeus et al. (2009), we approximate the turbulent energy cascade rate as

[TABLE]

where $D_{\rm K}=0.14$ is the cascade constant which was chosen to match the Kolmogorov constant. To study the effect of anisotropic turbulent cascade we simulated two regimes: one assumes the turbulent energy cascade is given by Eq. (9), the other assumes no cascade. The unperturbed spectrum of turbulence entering the free escape boundary at $z_{\mathrm{FEB}}$ is taken to be Kolmogorov, typically assumed for the interstellar medium. The incoming spectrum is normalized by

[TABLE]

3.2. Particle propagation model

With normal (i.e., non-Lévy-walk) diffusion the Monte Carlo simulation moves particles with a pitch-angle-scattering scheme that has been described in (Ellison et al., 1996). Briefly, after a time $\delta t$ much less than a gyroperiod a particle scatters isotropically and elastically in the local plasma frame through a angle $\delta\theta\leq\delta\theta_{\mathrm{max}}$, where $\delta\theta$ is chosen randomly between 0 and $\delta\theta_{\mathrm{max}}$. The maximum scattering angle is given by

[TABLE]

where $t_{c}=\lambda_{0}/v_{\mathrm{pf}}$ is the collision time, $\lambda_{0}(z,p)$ is the position and momentum dependent scattering length, $v_{\mathrm{pf}}$ is the particle speed in the rest frame of scattering centers, $r_{g}=pc/[eB_{\mathrm{ls}}(z,p)$] is the particle gyro-radius, and $B_{\mathrm{ls}}(z,p)$ is the local amplified magnetic field determined by summing fluctuations with wavelengths larger than $r_{g}$ (see Eq. (19) in reference Bykov et al., 2014, for a full description). In the normal scattering region we assume Bohm diffusion, i.e., $\lambda_{0}(z,p)=r_{g}(z,p)$.

For non-Lévy-walk scattering, particles always move a fraction of $\lambda_{0}(z,p)$ in the time interval $\delta t$, where $\lambda_{0}(z,p)$ is the mean free path a EP obtains scattering in the self-generated magnetic turbulence. We use the same Monte Carlo model as described in earlier work (see e.g., Section 2.7 of (Bykov et al., 2014)) to calculate MFA and $\lambda_{0}(z,p)$ except we now include the super-diffusion-induced mirror instability along with the resonant and non-resonant instabilities.

3.3. Super-diffusive EP propagation:

the Lévy-walk model

To model super-diffusive particle propagation we use a Lévy-walk model which assigns a random scattering length $\lambda_{\mathrm{LF}}$ to determine the path length. The probability density function of the $\lambda_{\mathrm{LF}}$ values has a power-law asymptotic form:

[TABLE]

This density function produces so-called “heavy tails” for $\nu\leq 3$, where $\nu=2$ is the well known Cauchy distribution. We describe below the specific algorithms to generate the random scattering lengths $\lambda_{\mathrm{LF}}$ with the power-law probability density functions. For completeness we discuss the Cauchy distribution in § 3.3.1 noting that since we only consider $\lambda_{\mathrm{LF}}>0$, $\nu=2$ in Eq. (12) gives the half-Cauchy distribution. We further note that while Cauchy distributions are used for mathematical convenience in semi-analytic calculations, they have an infinite mean and variance. The flexibility of the Monte Carlo model (and the fact that particles always have a speed less than $c$) allows us to derive a more general expression for the scattering length in § 3.3.2 for $2<\nu\leq 3$. We show examples with $\nu=2.1$ because this value gives a finite mean yet produces results similar to well-studied Cauchy ones.

3.3.1 The half-Cauchy distribution

With $\nu=2$ we assume the EP scattering length in the super-diffusive region $z_{\mathrm{FEB}}<z<z_{\mathrm{LF}}$ is determined by

[TABLE]

Here $\xi_{0}$ are random numbers uniformly distributed over the interval [0,1), where the brackets indicate values up to but not including 1. Values of $\xi_{0}\sim 1$ give extremely long scattering lengths, forcing a modification of Eq. (11), while $\lambda_{\mathrm{LF}}\sim 0$ for $\xi_{0}\sim 0$.

Integrating Eq. (13) from $p$ to $\infty$ we find the cumulative distribution function

[TABLE]

and the probability density function corresponding to Eq. (14) is

[TABLE]

3.3.2 Lévy-type power-law distributions

with $2<\nu\leq 3$

For $2<\nu\leq 3$ the recipe is somewhat more complicated. The scattering length in the Lévy-walk region of the shock precursor where $z<z_{\mathrm{LF}}$, is given by

[TABLE]

where $D_{0}=(\nu-1)/\nu$, and $\xi_{1}$ are random numbers uniformly distributed over the interval [0,1). The corresponding probability density function for $\lambda_{\rm LF}^{\ast}\equiv\lambda_{\rm LF}-\lambda_{0}$ is

[TABLE]

where

[TABLE]

(see Trotta & Zimbardo, 2015, for a full discussion).

The normalization and mean value are determined by

[TABLE]

and

[TABLE]

respectively. While Eq. (16) applies for $2<\nu\leq 3$, we restrict our calculations here to $\nu=2.1$ to ensure Eq. (20) yields a finite mean.

We implement Lévy-walk transport in a piecewise continuous way. In the precursor, downstream from some precursor position $z_{\mathrm{LF}}$, particles propagate diffusively with a mean free path $\lambda_{\mathrm{LF}}(z,p)=\lambda_{0}(z,p)$. Far upstream, between $z_{\mathrm{LF}}$ and the FEB at $z_{\mathrm{FEB}}$, we assume EPs propagate super-diffusively. Typically, $z_{\mathrm{LF}}\leq-10^{4}r_{g0}$, where $r_{g0}=m_{p}u_{0}c/(eB_{0})$ is the gyroradius of a proton with speed equal to the shock speed $u_{0}$ in the far upstream magnetic field $B_{0}$. For $z<z_{\mathrm{LF}}$, a particle is given a random scattering length $\lambda_{\mathrm{LF}}(z,p)$ generated with the Lévy-type stable distribution described below.

We restrict super-diffusive propagation to regions well upstream from the subshock where the self-generated turbulence has not reached saturation levels. Close to the FEB the turbulence is growing rapidly from the escaping EP flux and being convected downstream. Closer to the subshock (i.e., for $z>z_{\mathrm{LF}}$) the turbulence will be intense enough to destroy any long-range correlations and Bohm-like diffusion is assumed to occur. For $z<z_{\mathrm{LF}}$, we have a mixture of ballistic motion and diffusion, as described below.

3.4. Particle propagation with Lévy-walk

The Monte Carlo method we employ numerically solves a Boltzmann equation with a collision operator which is determined by collision frequencies averaged over the self-generated background turbulence. By replacing a “diffusion approximation” with a collision operator we are able to model pitch-angle-scattering controlled by short-scale fluctuations on the order of the particle gyroradius with arbitrary pitch angle distributions. Therefore, particle transport is not restricted to standard diffusion and super-diffusive, i.e., Lévy-walk, propagation can be directly modeled. All that is required for super-diffusion to be accurately modeled is to define a proper particle path length probability distribution.

If path lengths at a given position are determined by macroscopic, long-range field correlations on scales larger than the EP particle gyroradius, or by highly intermittent turbulence, then super-diffusion will occur. In the case when the path length distribution is Gaussian, or has a power-law index $\nu>3$, the probability of a EP having a free path that departs widely from the mean value is small and normal diffusion occurs with a well defined mean free path and mean square displacement $\left<\Delta z^{2}\right>=A_{D}t$.

On the other hand, if 2 $\leq\nu<3$ super-diffusive propagation occurs with $\left<\Delta z^{2}\right>=A_{D}t^{b}$ where $b>1$. In this case there is a non-negligible probability for the free path to be much longer than the mean. Physically, such a situation can be expected near the FEB where strong, unsaturated EP driven turbulence growth occurs. In this case, the turbulence is expected to be intermittent and long-range correlations are not immediately destroyed. Closer to the subshock, nonlinear interactions of strong magnetic fluctuations are likely to smooth out the intermittency and the downstream turbulence is likely to be statistically homogeneous. As a first approximation for this complicated situation, we model Lévy-walk propagation in the precursor in a region away from the subshock between $z_{\mathrm{FEB}}$ and $z_{\mathrm{LF}}$. The effect of varying $z_{\mathrm{LF}}$ is examined.

Once a path length distribution is specified, as with Eq. (3.3.2), the Monte Carlo algorithm determines $\left<\Delta z^{2}\right>(t)$ without further assumptions. As an illustration, we show in Fig. 1 Monte Carlo calculations with our algorithm of $\left<\Delta z^{2}\right>(t)$ vs. $t$ for two values of $\nu$, as indicated. This calculation is done in 3D geometry and the projection onto one axis is plotted. The super-diffusive case with $\nu=2.1$ yields a slope $b\simeq 1.76$. We note that in this example we restrict $\lambda_{\mathrm{LF}}$ to be equal to or greater than the Bohm limit, i.e., to the interval $[\lambda_{0},\infty)$. To test our Monte Carlo algorithm we performed the simulation without this restriction, where the scattering length is allowed to populate the interval $[0,\infty)$, and we obtained $b=1.86$ in good agreement with the scaling $b=4-\nu=1.9$ presented by (Zumofen & Klafter, 1993) and (Trotta & Zimbardo, 2015). The solid curve labeled $\nu=4.5$ is almost identical to the standard diffusion result shown with a dotted curve. For comparison we also show the ballistic case with $b=2$.

Particle transport with Lévy-walk is done in the following way. When a particle is at a position upstream from $z_{\mathrm{LF}}$ a random number $\xi_{1}$ is chosen and the particle’s scattering length, $\lambda_{\mathrm{LF}}(z,p)$, is found from Eq. (16). This determines the collision time $t_{c}=\lambda_{\mathrm{LF}}/v_{\mathrm{pf}}$. However, since $\lambda_{\mathrm{LF}}$ depends on $\lambda_{0}(p,z)$ and is position dependent in the nonlinear model, its value can change during $t_{c}$. To accommodate this we set the time interval

[TABLE]

during which the particle moves with a constant speed. After $\delta t(z,p)$, the particle scatters with a new $\lambda_{0}(z,p)$. The value of $\xi_{1}$ is kept fixed for $N_{\mathrm{scat}}$ scatterings where $N_{\mathrm{scat}}$ is 1 plus the integer part of $(\pi/A_{\theta})^{2}$. If a particle completes $N_{\mathrm{scat}}$ scatterings without escaping, the cycle is repeated with a new $\xi_{1}$. In the simulations presented below we set the parameter $A_{\theta}^{2}=\pi^{2}/2$. This large-angle scattering value was justified in (Ellison et al., 1996).

4. Results

In most of the following examples we calculate magnetic field amplification with super-diffusive EP propagation in nonlinear shocks where the energy and momentum conserving shock structure has been determined self-consistently. However, in order to isolate and highlight the effects of super-diffusion, we discuss some unmodified (UM) shocks with a discontinuous bulk-flow-velocity profile where energy and momentum are not conserved. For these unmodified shocks, $\lambda_{\mathrm{LF}}(z,p)=\lambda_{0}(p)$, i.e., spatially independent Bohm diffusion. These unmodified solutions are not, of course, intended to represent physical models.

In all cases, the shock speed $u_{0}=5000$ km s*-1*, the far upstream plasma density $n_{0}=0.3$ cm*-3*, the background magnetic field $B_{0}=3$ $\mu$G, and we accelerate only protons. The Fermi acceleration is limited by an upstream FEB at $z_{\mathrm{FEB}}=-10^{8}\,r_{g0}$, where $r_{g0}=m_{p}u_{0}c/(eB_{0})\simeq 5.6\times 10^{-9}$ pc.

As described in § 3.3, EPs move super-diffusively in the shock precursor between $z_{\mathrm{FEB}}$ and $z_{\mathrm{LF}}$ with a scattering length $\lambda_{\mathrm{LF}}(z,p)$ given by the Lévy-walk probability distribution Eq. (16) with $\nu=2.1$. The specific value $\nu=2.1$ is chosen because it is close to the Cauchy distribution and represents all of the features typical for prominent super-diffusive propagation while having a finite scattering length.

We show examples with $z_{\mathrm{LF}}=-10^{4}$, $-10^{5}$, and $-10^{6}\,r_{g0}$. For $z>z_{\mathrm{LF}}$, including downstream from the subshock, the EP scattering length is diffusive, i.e., $\lambda_{\mathrm{LF}}(z,p)=\lambda_{0}(z,p)$. For all nonlinear calculations, $\lambda_{0}(z,p)$ is determined with MFA from Bell’s instability, the resonant streaming instability, the non-resonant long-wavelength instability, and the mirror instability described here for the first time. For the parameters used here, the differences in $\lambda_{0}(z,p)$ derived with the additional mirror instability are modest except at the highest particle energies.

4.1. Nonlinear shock structure

In the top panels of Fig. 2 we show the nonlinear shock structure (solid black curve) in terms of the bulk plasma flow speed $u(x)$. The solid black curves in the middle and bottom panels show the momentum and energy fluxes for the self-consistent shock. For illustration, the dashed red curves show the corresponding quantities for an unmodified shock. An important element of the Monte Carlo simulation is that it contains a direct model of thermal leakage injection. The scattering assumptions described in § 2 determine the fraction of shock-heated particles that are injected into the Fermi acceleration mechanism. This, in turn, influences the overall acceleration efficiency in an internally self-consistent fashion.

As is clear from the dashed curves in Fig. 2, the Monte Carlo injection model predicts efficiencies that do not conserve energy and momentum in unmodified shocks. A consistent solution can be found without modifying the injection model by modifying the shock structure, as shown with the solid black curves in Fig. 2. As mentioned in § 3, the nonlinear bulk flow speed $u(z)$ is determined by iteration and results in momentum and energy conservation to within a few percent.

In this case, momentum and energy are conserved while still having a large Fermi acceleration efficiency. We emphasize that regardless of the injection process, shock modification must occur if Fermi acceleration is efficient. Furthermore, there must be a corresponding increase in the overall shock compression ratio, $R_{\mathrm{tot}}$, above the Rankine-Hugoniot value of $R_{\mathrm{tot}}\simeq 4$ for high Mach number shocks. The compression ratio is determined by the ratio of specific heats and the energy flux leaving the shock at the FEB. The distribution of escaping EPs is

[TABLE]

where $J(z_{\mathrm{FEB}},p)$ is the EP current at $z_{\mathrm{FEB}}$ measured in the upstream rest frame. Using the full anisotropy information provided by the Monte Carlo model, we define the position and momentum dependent EP current as

[TABLE]

where $f_{\mathrm{pf}}^{\mathrm{ep}}(z,p,\mu)$ is the distribution function of accelerated particles, per $d\mu$, in the local rest frame of the background plasma, $\mu=\cos\theta$, and $\theta$ is the angle between a particle’s momentum and the $z$-axis. The bottom panels in Fig. 2 show that $\sim 20\%$ of the far upstream energy flux $q_{\mathrm{esc}}=\int J(z_{\mathrm{FEB}},p)E(p)p^{2}dp$ (where $E(p)$ is the particle energy) escapes at $z=z_{\mathrm{FEB}}$ and the plot for $u(x)$ shows $R_{\mathrm{tot}}\simeq 7.2$. While $R_{\mathrm{tot}}$ increases above the test-particle Rankine-Hugoniot value, the subshock (indicated in the top right-hand panel) must decrease below the test-particle value. For the nonlinear shock in Fig. 2, $R_{\mathrm{sub}}\simeq 3$. These modifications to the shock structure from efficient diffusive shock acceleration produce the non-power-law behavior in the nonlinear distribution functions we discuss next.

4.2. Particle spectra

In Fig. 3 we show particle spectra measured downstream (DS) from the subshock and at the upstream FEB as indicated. The top panel shows spectra for unmodified shocks while spectra in the bottom panel are for self-consistent, nonlinear shocks, both with the super-diffusion parameter $\nu=2.1$ in Eq. (16). The proton distributions are calculated for $z_{\mathrm{LF}}=-10^{4}\,r_{g0}$ (black, solid and dashed curves) and without Lévy-flight transport (red, dot-dashed and blue, dotted curves). Super-diffusion is eliminated by placing $z_{\mathrm{LF}}$ outside of the FEB. In the unmodified examples, $\lambda_{0}$ was assumed to be spatially independent Bohm diffusion. The examples in Fig. 3 are calculated with cascading while those in Fig. 4 are calculated without cascading. We note that the statistical errors in Figs. 3 and 4 are small. Except for the escaping particles, variations are typically less than the line thickness.

These examples show that super-diffusion produces a high-energy cutoff in $f(z,p)$ that is broader than that for the diffusive case and occurs at a lower momentum. This is clearly reflected in the escaping distributions where the black dashed curves are with super-diffusion and the dotted blue curves are without. The broad cutoff results as particles with long Lévy-walk scattering lengths leave the FEB from deep within the precursor. It is significant that the broadening, while still present, is less in the self-consistent shocks. While not shown for clarity, we find that the cutoff is only weakly dependent on $z_{\mathrm{LF}}$.

4.3. Particle anisotropy with super-diffusion

The energetic particle current anisotropy, $A_{1}(z)$, is given by

[TABLE]

The partial anisotropy, $A_{2}(z,p)$, associated with the second spherical harmonic of the particle distribution is defined as

[TABLE]

and integrating over momentum gives

[TABLE]

where

[TABLE]

The anisotropies are both defined in the local plasma frame at all $z$.

In Fig. 5 we show dimensionless $A_{1}(z)$ and $A_{2}(z)$ for the unmodified shocks (top panel) and nonlinear shocks (bottom panel) discussed in Fig. 3. In the unmodified case, the bulk velocity profile is fixed with $R_{\mathrm{tot}}=4$ (see the dashed curves in Fig. 2) and the magnetic fluctuations spectra are also fixed, are position independent, and assume a Bohm-type scattering length $\lambda_{0}\propto p$. The nonlinear shocks are fully consistent in shock structure, $R_{\mathrm{tot}}$, self-generated magnetic turbulence, scattering length determination, and scattering center speed relative to the bulk speed of the background plasma (see § 4.5).

The important result here is that super-diffusive propagation with a Lévy-type distribution of particle scattering lengths results in second harmonics much stronger than produced with diffusive propagation. The magnitude of $A_{2}(z)$ is within a factor of a few to that of $A_{1}(z)$ in regions where super-diffusion is acting, while it is orders of magnitude less in regions with only diffusive propagation.

The second harmonic $A_{2}(z)$ is negative showing that the magnetic field partial pressure transverse to the local mean field is greater than the parallel partial pressure. This EP anisotropy will drive the so-called mirror instability, as we describe in more detail in Appendix A. The mirror instability is non-resonant where the growing magnetic fluctuations are nearly transverse wavevectors $2k_{\perp}^{2}>k_{\parallel}^{2}$ of scales larger than the gyroradius of particles, $R_{\mathrm{gp}}$, which are contributing into the transverse pressure, i.e., $k_{\perp}R_{\mathrm{gp}}<1$ but $k_{\perp}\lambda(p)>1$. The mirror mode is compressive and long-wavelength in the sense that $k_{\perp}R_{\mathrm{gp}}<1$. This characteristic may help increase the efficiency of scattering at the highest energy end of the accelerated particle spectrum.

We note that even though we emphasize effects produced for extremely efficient diffusive shock acceleration producing hard, concave spectra, as shown in Figs. 3 and 4, the growth rates we derive can be applied to weaker shocks. The mirror instability growth rate formulation can be applied to steep spectra and non-power-laws.

4.4. Magnetic field amplification

The energetic particle current and the quadrupole anisotropy, demonstrated in Fig. 5, drive resonant, short–, and long–wavelength instabilities, as well as the mirror instability first discussed here. This MFA is included in our nonlinear  model. The turbulence growth rates and the transport equations used were discussed in detail in sections 2.4 and 2.5 of (Bykov et al., 2014), while the growth rate for the mirror instability which is associated with the quadrupole anisotropy of EPs produced by super-diffusive propagation is presented in Appendix A.

The strong nonlinear aspects of MFA force a self-consistent description of the energy exchange between the EPs and the magnetic field, as well as with the bulk shock flow. We include the cascade of turbulence energy to shorter wavelengths taking into account the energy dissipation and heating of the background plasma. For comparison, we include models without turbulence cascade and show that cascading influences the magnetic fluctuation spectra, the total shock compression, and the downstream proton temperature.

In Fig. 6 we show the self-generated magnetic turbulence spectra, with and without cascading, at various positions relative to the subshock at $z=0$. These models include super-diffusion  beginning at $z_{\mathrm{LF}}=-10^{4}\,r_{g0}$. The top panel, without cascade, shows a strong spike in wave power at long wavelengths resulting from super-diffusion. With cascade, this turbulent energy is effectively shifted to shorter wavelengths.

4.5. Effective Scattering Center Velocity

The magnetic fluctuations produced by EP-driven instabilities in the shock precursor move relative to the bulk plasma with a speed $v_{\mathrm{scat}}(z)$. This is a highly nonlinear effect since $v_{\mathrm{scat}}(z)$ directly influences the effective compression ratio for diffusive shock acceleration which, in turn, determines the Fermi acceleration efficiency and the MFA. While virtually all work on diffusive shock acceleration attempting to consider a finite $v_{\mathrm{scat}}(z)$ assume it is some function of the Alfvén speed, $v_{\mathrm{alf}}$, there is no justification for this assumption other than ease of computation. A linear analysis of the four instabilities included in our model shows a vast variety of phase velocities. Significantly, even though the fastest growing EP-driven modes are highly anisotropic, the non-resonant Bell mode, as well as the long-wavelength and mirror modes, have phase speeds that are typically well below the local Alfvén speed calculated with the local amplified large-scale magnetic field (see Fig. 18 in Bykov et al., 2014).

The flexibility of the Monte Carlo method allows us to determine $v_{\mathrm{scat}}(z)$ from macroscopic energy conservation without making any assumptions regarding the amplified field. In Fig. 7 (bottom panels) we show $v_{\mathrm{scat}}(z)$ for the nonlinear cascade and non-cascade cases shown in Fig. 6.

As is clear from Fig. 7, $v_{\mathrm{scat}}$ is different in magnitude and spatial structure from $v_{\mathrm{alf}}$. Near the FEB $v_{\mathrm{scat}}>v_{\mathrm{alf}}$ but closer to the subshock $v_{\mathrm{scat}}$ can be orders of magnitude smaller. It is important to note that even though $v_{\mathrm{scat}}(z)$ may be small relative to both the bulk flow speed, it has a strong effect on the energy exchange between the accelerated particles, the bulk shock flow, and the magnetic field amplification. As in (Bykov et al., 2014), $v_{\mathrm{scat}}(z)$ is determined consistently with the shock structure modified by energetic particles and the magnetic field amplification.

5. Summary and conclusions

In collisionless astrophysical plasmas, particle transport is determined by charged particles interacting with magnetic turbulence and coherent magnetic structures over a broad wavelength range. These interactions are essential for collisionless shocks to form and accelerate particles to high energies. Details of the wave-particle interactions will influence energetic particles observed at Earth as well as radiative signatures of specific objects. While normal diffusion and advection have been proven to determine the long-range particle transport in quasi-homogeneous magnetic turbulence, studies both in the laboratory and in space plasmas have revealed a variety of sub- and super-diffusive regimes where particle transport may significantly depart from standard diffusive propagation for intermittent turbulence (e.g., during the growth of long-wavelength magnetic fluctuations in shock precursors).

We have presented the first nonlinear calculation of efficient Fermi shock acceleration that includes super-diffusion in a consistent manner. The Monte Carlo simulation we use, since it does not make a diffusion approximation, can model highly anisotropic particle distributions and magnetic field amplification and is well suited for these calculations. It includes nonlinear effects from thermal particle injection, shock modification, the self-generation of magnetic turbulence, turbulence cascade, and a consistently determined scattering center speed. While we fully expect that future large-scale PIC simulations will necessitate a refinement of our assumptions, accounting for important multi-dimensional effects is well beyond current PIC capabilities.

We show that super-diffusive particle transport in the shock precursor produces specific, anisotropic EP distributions which are characterized by a pronounced quadrupole anisotropy where the transverse particle pressure dominates the parallel particle pressure. We show that this type of anisotropy results in a mirror-like instability which is most prominent when the shock accelerated spectrum is $f(p)\propto p^{-4}$ or harder near the maximum energy of the accelerated particles (before the spectral break), as is expected for high Mach number shocks. Magnetic field and plasma structures produced by the mirror instability in cosmic plasmas were observed in planetary magnetosheaths (e.g., Hasegawa & Tsurutani (2011)) and may appear as well in high resolution Chandra images of young supernova remnants where synchrotron structures tracing the magnetic fields are prominent Tananbaum et al. (2014).

Super-diffusion also results in a broadening of the spectrum of EPs that escape from the shock precursor since there is a significant probability that EPs can leave the shock from deep inside the precursor. This broadening may impact models of $\gamma$-ray production by shocks interacting with dense molecular clouds and modify predictions for the integrated spectra of energetic particles.

ACKNOWLEDGMENTS A.M.Bykov and S.M.Osipov acknowledge support from RSF grant 16-12-10225. D.C.Ellison acknowledges support from NASA grant NNX11AE03G. Discussions concerning this work were held at the Aspen Center for Physics which is supported by National Science Foundation grant PHY-1066293.

Appendix A Mirror instability of anisotropic EP distributions

As mentioned in § 4.3, the Monte Carlo method follows an arbitrary angular distribution of accelerated particles and thus can determine the second spherical harmonic, $A_{2}(z,p)$, of the EP distribution at all positions across the nonlinear shock structure. A non-zero $A_{2}(z,p)$ indicates that there is a pressure anisotropy and this may result in the growth of a firehose or mirror instability depending on the ratio of partial pressures along and transverse to the mean magnetic field (e.g., Treumann & Baumjohann ????).

It has long been suggested that the mirror instability occurs in Maxwellian plasmas when the transverse pressure dominates the parallel pressure, as discussed by (Chandrasekhar et al. 1958; Sagdeev & Shafranov 1961; Hasegawa 1969; Southwood & Kivelson 1993; Hasegawa & Tsurutani 2011), and many others. Observational signatures of the mirror instability have been found in cosmic plasmas (e.g., Tsurutani et al. 2011) as well. However, in highly nonthermal systems like the collisionless shocks discussed in this paper the non-thermal relativistic particle pressure is important. In the downstream flow it is comparable to the thermal pressure and it can be well above the thermal particle pressure in the shock precursor.

Therefore, consistent nonlinear shock solutions must determine the growth rate of the mirror instability, in the thermal background plasma, driven by a superthermal particle pressure anisotropy that occurs with super-diffusive  transport. Our derivation treats the accelerated particles, the source of free energy for the mirror instability, kinetically, while the background plasma is described with an MHD approximation (e.g., Bykov et al. 2013). Since injection is a continuous process in the Monte Carlo simulation, all particles that make at least one crossing of the subshock (at $z=0$) from downstream to upstream are included as superthermal.

For the background plasma we have

[TABLE]

where $\rho$, $\mathbf{u}$, and $p_{g}$ are the background plasma density, macroscopic velocity, and pressure respectively. In addition, the continuity equation is

[TABLE]

the electric and magnetic fields $\mathbf{E}$ and $\mathbf{B}$ satisfy

[TABLE]

and

[TABLE]

where $f^{\mathrm{ep}}$, $n_{\mathrm{ep}}$, and $\mathbf{j}^{\mathrm{ep}}$ are the EP distribution function, EP concentration, and electric current, $\mathcal{E}$ is the EP particle energy, $\displaystyle\widehat{\mathbf{\mathcal{O}}}=\mathbf{p}\times\frac{\partial}{\partial\mathbf{p}}$ is the momentum rotation operator, $c$ is the speed of light, $e$ is the particle charge, and $\displaystyle\Omega=eB_{0}c/\mathcal{E}$ is the EP particle gyro-frequency.

The right-hand-side of Eq. (A4), $I[f^{\mathrm{ep}}]$, is the collision operator describing EP interactions with magnetic fluctuations carried by the background plasma. In the background plasma frame,

[TABLE]

where $f^{\mathrm{ep}}_{\mathrm{iso}}$ is the isotropic part of the distribution function, $\nu_{c}=\epsilon\Omega$ is the EP scattering frequency by magnetic fluctuations with $\epsilon\leq 1$, and $\int\mathbf{p}I[f^{\mathrm{ep}}]d^{3}p=-\epsilon B_{0}\mathbf{j}^{\mathrm{ep}}/c$ in Eq. (A).

In the rest frame of the background plasma we impose small perturbations of the local plasma parameters in Eq. (A)–(A4) as $\xi=\xi_{0}+\delta\xi$, with $\delta\xi\propto\exp(i\mathbf{k}\mathbf{r}-i\omega t)$. Assuming the adiabatic equation for the background plasma, $\nabla\delta p_{g}=a_{0}^{2}\nabla\delta\rho$, with $a_{0}=\sqrt{\gamma_{g}p_{0}/\rho_{0}}$, where $\gamma_{g}$ is the adiabatic index. Furthermore, we define $\mathbf{e}_{z}=\mathbf{B}_{0}/B_{0}$ as the parallel direction, $\mathbf{e}_{x}$ is the transverse direction, and $\mathbf{k}=k_{\parallel}\mathbf{e}_{z}+k_{\perp}\mathbf{e}_{x}$.

Since we are only concerned here with instabilities due to the anisotropic EP pressure (i.e., the quadrupole anisotropy of the EP distribution), we consider the unperturbed EP distribution function with no mean electric current (see Schure et al. 2012; Bykov et al. 2013; Marcowith et al. 2016, for discussions of the resonant and current driven instabilities). With no mean current, the quadrupole anisotropy has the form

[TABLE]

where $\int_{0}^{\infty}N(p)p^{2}dp=1$ and $\chi$ is a quadrupole anisotropy parameter with $|\chi|<1$.

Keeping only linear responses to the perturbations in Eq. (A)–(A2), we obtain

[TABLE]

and

[TABLE]

where $v_{a}=B_{0}/\sqrt{4\pi\rho_{0}}$ is the Alfvén velocity. The linearized Eq. (A4) has the form

[TABLE]

where $\theta$ and $\varphi$ are the pitch and azimuthal angles between the EP particle velocity and the direction $\mathbf{e}_{z}$, correspondingly.

We first consider the weakly collisional case with $\epsilon\ll 1$. In the long-wavelength regime, $k_{\parallel}v/\Omega\ll 1$, $k_{\perp}v/\Omega\ll 1$, and $\omega/\Omega\ll 1$. Then one obtains the response of the superthermal particle current $\delta\mathbf{j^{ep}}$ on the magnetic field perturbation $\delta\mathbf{B}$ in the form

[TABLE]

and

[TABLE]

Consider the brace in Eq. (A). With the assumption for the unperturbed quadrupole anisotropy given by Eq. (A6), the second and third terms in the brace are of the same order of magnitude. The second term gives the well known firehose instability which grows if $\chi>0$, while the third term results in the mirror instability if $\chi<0$. Using Eqs. (A6) and (25), the parameter $\chi$ can be connected to $A_{2}(p)$ with $A_{2}(p)=n_{\mathrm{ep}}N(p)\chi/5$. Note that if the EP scattering rate by magnetic turbulence $\nu_{c}\sim\Omega$, which may occur if $\epsilon\sim$ 1, then the mirror instability is suppressed.

Integrating Eq. (A) over $\mu$, with account taken of Eq. (A6), one obtains

[TABLE]

where $\tau=k_{\parallel}v/\omega$ and $\Theta(z)$ is the Heaviside step function. The asymptotic form of Eq. (A) for $\tau\ll 1$ is

[TABLE]

while, for $\tau\gg 1$,

[TABLE]

The response of the electric current carried by the energetic particles in these limits is

[TABLE]

and

[TABLE]

Let us define the unperturbed pressure

[TABLE]

and note that for the distribution function Eq. (A6),

[TABLE]

[TABLE]

and

[TABLE]

Now, substituting Eqs. (A15) and (A16) into (A) and (A8) we obtain the dispersion relations

[TABLE]

and

[TABLE]

The dispersion equation splits into two independent equations where the mirror instability is determined by Eq. (A) under conditions $2k_{\perp}^{2}>k_{\parallel}^{2}$ and $\chi<0$. For $k_{\perp}^{2}\gg k_{\parallel}^{2}$ Eq. (A) simplifies to

[TABLE]

Here the perturbation of the magnetic field is mostly $\delta B_{z}$ which are connected to the background plasma density variations by $\displaystyle\delta\rho\approx\rho_{0}\frac{\delta B_{z}}{B_{0}}$.

Then from Eq. (A20), with account taken that $k\approx k_{\perp}$, one finally obtains

[TABLE]

Eq. (A24) shows that growing modes occur when $\left(v_{a}^{2}+a_{0}^{2}+2\frac{\delta P}{\rho_{0}}\right)<0$, which may happen only if $\delta P\propto\chi<0$, i.e., the mirror instability is driven by anisotropic EP pressure.