Stochastic Ion Acceleration by the Ion-cyclotron Instability in a Growing Magnetic Field

TL;DR

This study demonstrates that ions can be stochastically accelerated via the ion-cyclotron instability in growing magnetic fields, especially at low plasma beta, with implications for accretion flows around black holes.

Contribution

It reveals the conditions under which ion acceleration occurs through the ion-cyclotron instability in growing magnetic fields, highlighting the role of plasma beta and magnetic field amplification.

Findings

Ion acceleration produces a non-thermal tail in the ion energy spectrum.

Efficient acceleration requires initial plasma beta less than or approximately 1.

The process is applicable in sub-relativistic regimes relevant to black hole accretion flows.

Abstract

Using 1D and 2D particle-in-cell (PIC) simulations of a plasma with a growing magnetic field $\vec{B}$, we show that ions can be stochastically accelerated by the ion-cyclotron (IC) instability. As $\vec{B}$ grows, an ion pressure anisotropy $p_{\perp,i} > p_{||,i}$ arises, due to the adiabatic invariance of the ion magnetic moment ($p_{||,i}$ and $p_{\perp,i}$ are the ion pressures parallel and perpendicular to $\vec{B}$). When initially $\beta_i = 0.5$ ($\beta_i \equiv 8\pi p_i/|\vec{B}|^2$, where $p_i$ is the ion isotropic pressure), the pressure anisotropy is limited mainly by inelastic pitch-angle scattering provided by the IC instability, which in turn produces a non-thermal tail in the ion energy spectrum. After $\vec{B}$ is amplified by a factor $\sim 2.7$, this tail can be approximated as a power-law of index $\sim 3.4$ plus two non-thermal bumps, and accounts for $2-3\%$ of the ions and $\sim 18\%$ of their kinetic energy. On the contrary, when initially $\beta_i =2$, the ion scattering is dominated by the mirror instability and the acceleration is suppressed. This implies that efficient ion acceleration requires that initially $\beta_i \lesssim 1$. Although we focus on cases where $\vec{B}$ is amplified by plasma shear, we check that the acceleration occurs similarly if $\vec{B}$ grows due to plasma compression. Our results are valid in a sub-relativistic regime where the ion thermal energy is $\sim 10\%$ of the ion rest mass energy. This acceleration process can thus be relevant in the inner region of low-luminosity accretion flows around black holes.

Tables

Table 1. Table 1 Parameters of the simulations

Runs	${\omega }_{c,i}^{\text{init}}/s$	$m_i/m_e$	${\beta }_i^{\text{init}}$	N${}_{\text{ppc}}$	L/R${}_{}{}^{\text{init}}{}_{L,i}{}^{}$[D]
	(or ${\omega }_{c,i}^{\text{init}}/q$)
S2m2b0.5	800	2	0.5	160	60[2]
S2m2b2	800	2	2	160	60[2]
S2m10b0.5	800	10	0.5	160	60[2]
S2m10b2	800	10	2	160	60[2]
S1m2b0.5	800	2	0.5	640	150[1]
S1m2b2	800	2	2	640	150[1]
S1m8b0.5	800	8	0.5	640	150[1]
S1m32b0.5	800	32	0.5	640	150[1]
S1m128b0.5	800	128	0.5	1280	150[1]
S1m8b0.5b	400	8	0.5	640	150[1]
S1m8b0.5c	1600	8	0.5	640	150[1]
S1m8b0.5d	3200	8	0.5	3600	150[1]
C1m8b0.5a	1600	8	0.5	3600	190[1]
C1m8b0.5b	3200	8	0.5	3600	190[1]
C1m16b0.5a	1600	16	0.5	3600	190[1]
C1m16b0.5b	3200	16	0.5	3600	190[1]

Equations

$$\frac{d{U}_j}{dt}=r\Delta p_j,$$

$${\nu }_{eff,i}\propto \frac{{\omega }_{c,i}}{B_0^2{\gamma }_i}\frac{d(\delta B_z^2)}{d\ln (k)}.$$

$$\frac{\partial p_{\parallel ,i}}{\partial t}+\nabla \cdot (p_{\parallel ,i}\text{v})+2p_{\parallel ,i}\hat{\text{b}}\hat{\text{b}}:\nabla \text{v}=\frac{2}{3}{\nu }_{eff,i}\Delta p,$$

$${\nu }_{eff,i}\approx \frac{3}{2}s(p_{\parallel ,i}/\Delta p_i).$$

$$\frac{\omega }{k}-v_{||}=\frac{{\omega }_{c,i}}{{\gamma }_{i}k},$$

$${\gamma }_i(kc/{\omega }_{c,i})\approx 2.$$

$${\gamma }_i(kc/{\omega }_{c,i})\approx \mathrm{1.5.}$$

$${\nu }_{eff,i}\approx \frac{3}{4}q(p_{\parallel ,i}/\Delta p_i).$$

$$x^{\prime }=x,y^{\prime }=\Gamma (y-vt),z^{\prime }=z,\text{and}{t}^{\prime }=\Gamma (t-vy/c^2),$$

$$x_1^{\prime }\equiv x^{\prime }-\frac{y^{\prime }s{t}^{\prime }}{1+s^2{t}^{\prime 2}}.$$

$$\frac{\partial }{\partial x}=\frac{1}{1+s^2{t}^2}\frac{\partial }{\partial x_1^{\prime }}+st\frac{\partial }{\partial y^{\prime }},\frac{\partial }{\partial y}=\frac{-st}{1+s^2{t}^2}\frac{\partial }{\partial x_1^{\prime }}+\frac{\partial }{\partial y^{\prime }},\text{and}\frac{\partial }{\partial z}=\frac{\partial }{\partial z^{\prime }}.$$

$$\frac{\partial B_{x}^{\prime}}{\partial t^{\prime}}=c\frac{st^{\prime}}{1+s^{2}t^{\prime 2}}\frac{\partial E_{z}^{\prime}}{\partial x^{\prime}_{1}},\hskip 10.00002pt\frac{\partial B_{y}^{\prime}}{\partial t^{\prime}}=c\frac{1}{1+s^{2}t^{\prime 2}}\frac{\partial E_{z}^{\prime}}{\partial x^{\prime}_{1}}-sB_{x}^{\prime},\hskip 10.00002pt\frac{\partial B_{z}^{\prime}}{\partial t^{\prime}}=\frac{-c}{1+s^{2}t^{\prime 2}}\Big{(}\frac{\partial E_{y}^{\prime}}{\partial x^{\prime}_{1}}+st^{\prime}\frac{\partial E_{x}^{\prime}}{\partial x^{\prime}_{1}}\Big{)},$$

$$\frac{\partial E_{x}^{\prime}}{\partial t^{\prime}}=-c\frac{st^{\prime}}{1+s^{2}t^{\prime 2}}\frac{\partial B_{z}^{\prime}}{\partial x^{\prime}_{1}}-4\pi J^{\prime}_{x},\hskip 10.00002pt\frac{\partial E_{y}^{\prime}}{\partial t^{\prime}}=-c\frac{1}{1+s^{2}t^{\prime 2}}\frac{\partial B_{z}^{\prime}}{\partial x^{\prime}_{1}}-sE_{x}^{\prime}-4\pi J^{\prime}_{y},\hskip 10.00002pt\frac{\partial E_{z}^{\prime}}{\partial t^{\prime}}=\frac{c}{1+s^{2}t^{\prime 2}}\Big{(}\frac{\partial B_{y}^{\prime}}{\partial x^{\prime}_{1}}+st^{\prime}\frac{\partial B_{x}^{\prime}}{\partial x^{\prime}_{1}}\Big{)}-4\pi J^{\prime}_{z}.$$

$$\frac{d\textbf{{p}}^{\prime}}{dt^{\prime}}=q_{c}\big{(}\textbf{{E}}^{\prime}+\frac{\textbf{{u}}^{\prime}}{c}\times\textbf{{B}}^{\prime}\big{)}+sp_{x}^{\prime}\hat{y},$$

$$\frac{d{x}_1^{\prime }}{d{t}^{\prime }}=u_x^{\prime }\frac{1}{1+s^2{t}^{\prime 2}}-u_y^{\prime }\frac{s{t}^{\prime }}{1+s^2{t}^{\prime 2}},$$

$$\frac{\partial {\rho }_c^{\prime }}{\partial t^{\prime }}+\frac{\partial J_x^{\prime }}{\partial x^{\prime }}+\frac{\partial (J_y^{\prime }+s{t}^{\prime }{J}_x^{\prime })}{\partial y^{\prime }}+\frac{\partial J_z^{\prime }}{\partial z^{\prime }}=0,$$

$$\frac{d{x}^{\prime }}{d{t}^{\prime }}=u_x^{\prime },\frac{d{y}^{\prime }}{d{t}^{\prime }}=u_y^{\prime }+s{t}^{\prime }{u}_x^{\prime },\text{and}\frac{d{z}^{\prime }}{d{t}^{\prime }}=u_z^{\prime },$$

$$\frac{\partial\rho_{c}^{\prime}}{\partial t^{\prime}}+\frac{\partial}{\partial x^{\prime}_{1}}\Big{(}\frac{J_{x}^{\prime}}{1+s^{2}t^{\prime 2}}-\frac{st^{\prime}J_{y}^{\prime}}{1+s^{2}t^{\prime 2}}\Big{)}=0.$$

Full text

1. Introduction

Stochastic (or second-order Fermi) acceleration by plasma turbulence is considered a viable mechanism for producing non-thermal particles in several astrophysical environments. This process can in principle be driven by MHD plasma waves (e.g., Chandran, 2003; Cho & Lazarian, 2006; Lynn et al., 2014) and by kinetic plasma modes (e.g., Dermer et al., 1996; Petrosian & Liu, 2004). In this work we use particle-in-cell (PIC) plasma simulations to show that ions can be stochastically accelerated by ion-cyclotron (IC) waves driven unstable in the presence of an ion pressure anisotropy with $p_{\perp,i}>p_{||,i}$ (where $p_{\perp,i}$ and $p_{||,i}$ are the ion pressures perpendicular and parallel to the local magnetic field B, respectively).

The condition $p_{\perp,i}\neq p_{||,i}$ is naturally expected in turbulent, weakly collisional plasmas. In these environments, Coulomb collisions are not able to break the adiabatic invariance of the magnetic moment $\mu_{i}$ of ions, which is defined as $\mu_{i}\equiv v_{\perp,i}^{2}/B$, where $v_{\perp,i}$ is the ion velocity perpendicular to B and $B=|\textbf{{B}}|$. Thus, if $B$ grows (decreases), the conservation of $\mu_{i}$ will naturally produce a pressure anisotropy with $p_{\perp,i}>p_{||,i}$ ($p_{\perp,i}<p_{||,i}$). Examples of weakly collisional astrophysical plasmas where the condition $p_{\perp,i}\neq p_{\parallel,i}$ is possible are low-luminosity accretion disks around compact objects (e.g., Sharma et al., 2006), the intracluster medium (ICM; Schekochihin et al., 2005; Lyutikov, 2007), and the heliosphere (e.g., Bale et al., 2009; Maruca et al., 2011; Verscharen et al., 2019).

In these systems, the growth of ion pressure anisotropy is expected to be regulated by kinetic instabilities, which break the adiabatic invariance of $\mu_{i}$ via pitch-angle scattering of the ions. In the $p_{\perp,i}>p_{\parallel,i}$ regime, there are two relevant instabilities: the mirror and the ion-cyclotron (IC) instabilities. The mirror instability consists of non-propagating, compressional modes, with their dominant modes having wave vectors k oblique to the direction of B (Hasegawa, 1969; Southwood & Kivelson, 1993). The IC instability, on the other hand, consists of propagating electromagnetic modes, with their dominant waves having $\textbf{{k}}\parallel\textbf{{B}}$ (Anderson et al., 1991; Gary, 1992). Whether the ion pitch-angle scattering is dominated by the IC or mirror instability essentially depends on how fast the instabilities grow for a given plasma regime. In this work we show that, in a regime dominated by the IC instability, significant non-thermal ion acceleration can occur due to scattering by the IC waves.

Our study will use particle-in-cell (PIC) plasma simulations to study a homogeneous plasma in which ion and electron pressure anisotropies are self-consistently produced by the continuous growth of a background magnetic field $B$. In our simulations, the magnetic field will grow on time scales significantly longer than the initial, exponential growth regime of the mirror and IC instabilities. This will allow us to capture the long-term, saturated state of the instabilities, which should be the dominant regime in astrophysical systems where $B$ experiences significant amplifications. In most of our simulations, the magnetic field will be amplified through imposing a slow shear motion in the plasma, which will increase $B$ due to magnetic flux conservation. However, we will also use simulations of compressing plasmas to show that our main results are fairly independent of the specific mechanism that drives the growth of $B$. In this study we focus on conditions applicable to the inner regions of low-luminosity accretion flows around black holes. This will be done by assuming in all of our runs a hot plasma with initially equal ion and electron temperatures ($T_{i}=T_{e}$) and with $k_{B}T_{i}/m_{i}c^{2}=0.05$ ($k_{B}$ is the Boltzmann constant, $m_{i}$ is the mass of the ions, and $c$ is the speed of light).

Our paper is organized as follows. In §2 we present our simulation method and setup. In §3 we use 2D simulations to show that ions can be accelerated by the IC instability under the condition that initially $\beta_{i}\lesssim 1$ ($\beta_{i}\equiv 8\pi p_{i}/|\textbf{{B}}|^{2}$, where $p_{i}$ is the ion isotropic pressure). In §4 we use 1D simulations to clarify the role of the IC and mirror modes in the acceleration, as well as to show that our results are independent of the rate at which $B$ is amplified and of the numerical ion to electron mass ratio, $m_{i}/m_{e}$. In §5 we describe the acceleration in further details, connecting the growth of IC modes of different wavenumber $k$ with the acceleration of ions of different energy. In §6 we show that our results are fairly independent on whether $B$ is amplified via plasma shear or compression. In §7 we summarize our results and present our conclusions. Additionally, in Appendix §A we provide details on the implementation of our 1D simulations, and in Appendix §B we use linear theory to analyze the applicability of our simulation results to realistic astrophysical environments.

2. Simulation Setup

We use the PIC code TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005) to simulate both a shearing and a compressing plasma made of ions and electrons. In the shearing case, the plasma is initially in presence of a homogeneous initial magnetic field that points along the $x$ axis, $\textbf{{B}}=B_{0}\hat{x}$. This field is amplified by imposing a shear plasma velocity $\textbf{{v}}=-sx\hat{y}$ (represented by red arrows in Fig. 1$a$), where $x$ is the distance along $\hat{x}$ and $s$ is the shear rate. This way the background magnetic field B in the simulation permanently increases and changes direction due to magnetic flux conservation, with its $y$-component evolving as $dB_{y}/dt=-sB_{0}$, while $dB_{x}/dt=dB_{z}/dt=0$ (Fig. 1$b$ shows how B changes orientation for $t>0$). Due to $\mu_{j}$ conservation, this magnetic growth drives $p_{\perp,j}>p_{||,j}$ during the whole simulation, allowing the triggering of kinetic instabilities that limit the pressure anisotropies.

Our 2D shearing runs use initially square simulation domains (as the one depicted in Fig. 1$a$) that follow the mean shear motion of the plasma. Therefore, the 2D domain acquires a parallelogram shape for $t>0$ (as shown in Fig. 1$b$). The positions of the plasma particles are therefore given in terms of the so called ‘shearing coordinates’, which are described both in Appendix A of this paper as well as in the Appendix of Riquelme et al. (2012). In our 1D shearing runs, on the other hand, the simulation domain corresponds to the blue narrow stripe shown in Figs. 1$a$ and 1$b$, which also moves with the shearing flow. Since the symmetry axis of this 1D domain is permanently parallel to B, our 1D approach allows to capture waves that propagate parallel to B. The self-consistent implementation of the 1D runs requires a small change in the definition of our shear coordinates, which is explained in detail in Appendix A. In §4.1 we show that our 1D runs give essentially the same results as our 2D runs as long as the dominant instabilities produce modes parallel to B, such as the IC instability.

In our compressing plasma runs, on the other hand, the simulation box is compressed along the two directions perpendicular to the background field B, producing both the permanent growth of B and of $p_{\perp,j}/p_{\parallel,j}$. For this we use the same setup as in Sironi & Narayan (2015). In this case, B evolves as $\textbf{{B}}=\hat{x}B_{0}/(1+qt)^{2}$, where the constant $q$ provides the time scale for the plasma compression.

Our plasma parameters are the initial temperature of ions and electrons ($T_{i}$ and $T_{e}$), the initial ratio between ion pressure and magnetic pressure ($\beta_{i}^{\textrm{init}}$), the ion to electron mass ratio $m_{i}/m_{e}$, and the ion “magnetization”, which is defined as the ratio between the initial ion cyclotron frequency ($\omega_{c,i}^{\textrm{init}}$) and $s$ (for shearing plasma runs) or $q$ (for compressing plasma runs). The initial ion cyclotron frequency is defined as $\omega_{c,i}^{\textrm{init}}=eB_{0}/m_{i}c$, with $e$ and $B_{0}$ being the magnitude of the electron and ion electric charges and the initial magnetic field.

As mentioned above, all of our shearing and compressing simulations have initially $T_{i}=T_{e}$ and $k_{B}T_{i}/m_{i}c^{2}=0.05$. Also, these runs use $m_{i}/m_{e}$ and ion magnetizations much smaller than expected in real astrophysical settings.111 For instance, at $\sim 10$ Schwarzschild radii from the super-massive black hole Sgr A*, one expects $\omega_{c,i}/s\sim 10^{8}$ (e. g., Ponti et al., 2017), where we have approximated $s$ as the Keplerian angular frequency at that radius. Because of this, the dependence of the ion acceleration on these parameters will be carefully assessed. (With our main conclusion being that neither $m_{i}/m_{e}$ or $\omega_{c,i}^{\textrm{init}}/s$ (or $\omega_{c,i}^{\textrm{init}}/q$) play a significant role.)

The numerical parameters in our runs are: the number of macro-particles per cell (N${}_{\textrm{ppc}}$), the electron skin depth in terms of grid point spacing ($c/\omega_{p,e}/\Delta_{x}$, where $\omega_{p,e}^{2}=4\pi n_{e}e^{2}/m_{e}$ is the electron plasma frequency and $n_{e}$ is the electron number density), and the box size in terms of the initial ion Larmor radius ($L/R_{L,i}^{\textrm{{init}}}$; $R_{L,i}^{\textrm{{init}}}=v_{th,i}/\omega_{c,i}^{\textrm{init}}$, where $v_{th,i}^{2}=k_{B}T_{i}/m_{i}$). Table 1 shows a summary of our key simulations. We ran a series of simulations ensuring that the numerical parameters do not affect our results. The runs used just for numerical convergence are not in Table 1.

3. Ion Acceleration by the IC Instability

We use 2D, shearing plasma simulations to show that ions can be stochastically accelerated by the IC modes. First, we show the example of runs dominated by the IC instability (with $\beta_{i}^{\textrm{init}}=0.5$), demonstrating that in this case a prominent non-thermal tail appears. Then, using simulations with $\beta_{i}^{\textrm{init}}=2$, we show that for $\beta_{i}^{\textrm{init}}\gtrsim 1$ the mirror instability dominates, with a corresponding suppression of the accelerating effect of the IC modes. We use 2D simulations with $m_{i}/m_{e}=2$ and 10 to show that, as long as the ion physics is concerned, our results are fairly independent of the value of $m_{i}/m_{e}$.

3.1. IC vs. mirror dominated regimes

Figure 2$a$ shows in blue-solid and red-solid lines the respective evolutions of $p_{\perp,i}$ and $p_{\parallel,i}$ for run S2m2b0.5, which uses $\beta_{i}^{\textrm{init}}=0.5$, $m_{i}/m_{e}=2$ and $\omega_{c,i}^{\textrm{init}}/s=800$. We see that until $t\cdot s\approx 1$ the evolutions of both $p_{\perp,i}$ and $p_{\parallel,i}$ are in agreement with the “double adiabatic” prediction (dashed lines; Chew et al., 1956), which is due to the conservation of $\mu_{i}$ and of the second adiabatic invariant. At $t\cdot s\gtrsim 1$, the adiabatic evolution of $p_{\perp,i}$ and $p_{\parallel,i}$ is broken by the appearance of ion pressure anisotropy instabilities, which produce rapid pitch-angle scattering of the ions.

We see from Figure 3 that these instabilities are dominated by IC modes. Indeed, Figs. 3$a$, 3$b$ and 3$c$ show a snapshot at $t\cdot s=2$ of the three components of $\delta\textbf{{B}}$ (where $\delta\textbf{{B}}\equiv\textbf{{B}}-\langle\textbf{{B}}\rangle$ and $\langle\rangle$ denotes an average over the entire box volume). Considering that the black arrows represent the direction of $\langle\textbf{{B}}\rangle$, we see that $\delta\textbf{{B}}$ is dominated by nearly parallel modes, which mainly appear in $\delta B_{x}$ and $\delta B_{z}$. This is indeed consistent with the presence of transverse, circularly polarized IC modes. $\delta B_{y}$ shows a mixture of the nearly parallel modes, plus subdominant oblique modes, which are consistent with the presence of mirror modes. Indeed, these modes mainly contribute to the $\delta\textbf{{B}}$ components parallel to the plane of the simulation (they show no $\delta B_{z}$ component), which is in line with the expectation that $\delta\textbf{{B}}$ of the mirror modes is nearly perpendicular to $\textbf{{k}}\times\textbf{{B}}$ (Pokhotelov et al., 2004).

The nearly parallel IC modes can also be seen from Figs. 3$d$, 3$e$ and 3$f$, which show the fluctuations in the electric field, $\delta\textbf{{E}}$.222Since our simulations are performed in the “shearing coordinate” frame (Riquelme et al., 2012), there is not an electric field associated to the large-scale shearing motion of the plasma. Therefore, $\delta\textbf{{E}}$ corresponds to the entire electric field present in the simulation. This electric field is expected since the IC modes have finite phase velocities, $v_{\phi}$, which is related to $\delta\textbf{{E}}$ by $|\delta\textbf{{E}}|=|\delta\textbf{{B}}|v_{\phi}/c$ (this is a consequence of Faraday’s law applied to the transverse IC modes). The mirror modes, on the other hand, are ‘purely growing’ (see, e.g., Southwood & Kivelson, 1993), which means that their phase velocity vanishes. This implies that, as we see in Fig. 3, no electric field associated to the subdominant mirror modes should be present.

The dominance of the IC modes can also be seen from Figure 4$a$, which shows in solid lines the magnetic energy of $\delta\textbf{{B}}$ along different axes as a function of time for run S2m2b0.5. This energy is expressed in terms of the $\delta\textbf{{B}}$ components parallel to $\langle\textbf{{B}}\rangle$ ($\delta B_{\parallel}$; solid-blue), perpendicular to $\langle\textbf{{B}}\rangle$ but parallel to the plane of the simulation ($\delta B_{xy,\perp}$; solid-red), and perpendicular to both $\langle\textbf{{B}}\rangle$ and the plane of the simulation ($\delta B_{z}$; solid-green). During most of the simulation, the energy of the magnetic fluctuations is indeed contained mainly in $\delta B_{xy,\perp}$ and $\delta B_{z}$, implying that the IC modes have the largest amplitude during most of the simulation time. By the end of the run ($t\cdot s=2.5$), however, $\delta B_{\parallel}^{2}$ becomes comparable to $\delta B_{xy,\perp}^{2}$ and $\delta B_{z}^{2}$, implying that in the long term the mirror fluctuations can still reach amplitudes comparable to the IC modes. We have thus decided to concentrate on the regime where the IC instability clearly dominates the pitch-angle scattering of the ions by running the simulations until $t\cdot s=2.5$ (thus with a maximum $B$ amplification factor of $\sim 2.7$).

Figure 4$b$, on the other hand, shows the evolution of the same magnetic energy components for the 2D run S2m2b2 ($m_{i}/m_{e}=2$ and $\beta_{i}^{\textrm{init}}=2$). In this run, the ions are under the same conditions as in run S2m2b2, but with a smaller initial background magnetic field so that $\beta_{i}^{\textrm{init}}=2$. We see that in this case $\delta B_{z}$ and $\delta B_{xy,\perp}$ are subdominant in the saturated stage of the instabilities, and the energy in the magnetic fluctuations is dominated by $\delta B_{\parallel}^{2}$. This result indicates that the oblique mirror modes are more prominent than the IC modes in this case, with the transition from IC-dominated to mirror-dominated regimes happening at $\beta_{i}^{\textrm{init}}\sim 1$.

In order to explore the sensitivity of this transition to $m_{i}/m_{e}$, in Figs. 4$a$ and 4$b$ we overplot $\delta B_{\parallel}^{2}$, $\delta B_{xy,\perp}$, and $\delta B_{z}$ for simulations S2m10b0.5 and S2m10b2, which have the same ion conditions as in runs S2m2b0.5 and S2m2b2 (i.e., the same values of $\beta_{i}^{\textrm{init}}$, $k_{B}T_{i}/m_{i}c^{2}$ and $\omega_{c,i}^{\textrm{init}}/s$), but with $m_{i}/m_{e}=10$ instead of $m_{i}/m_{e}=2$. We see that for the two $\beta_{i}^{\textrm{init}}$ the evolutions of $\delta B_{xy,\perp}$, $\delta B_{z}$, and $\delta B_{\parallel}$ are fairly independent of $m_{i}/m_{e}$. Thus the mass ratio does not appear to affect significantly the dominance of the IC instability for $\beta_{i}^{\textrm{init}}\lesssim 1$.

Additionally, in Appendix §B we use linear theory calculations to show that the condition $\beta_{i}^{\textrm{init}}\lesssim 1$ for the dominance of the IC instability should continue to hold even in realistic astrophysical plasma conditions, with $m_{i}/m_{e}=1836$, $\omega_{c,i}^{\textrm{init}}/s\gg 800$ and $T_{e}=T_{i}$. In the next section we show that this IC dominance also results in a significant non-thermal ion acceleration, which is strongly suppressed when the mirror modes dominate.

3.2. Ion heating and acceleration

It is well known that in a collisionless, shearing plasma the particles are heated by the so called “anisotropic viscosity”. Indeed, for a homogeneous plasma subject to shear, the internal energy density for species $j$, $U_{j}$, evolves as (Kulsrud, 1983; Snyder et al., 1997)

[TABLE]

where $r$ is the growth rate of the magnetic field ($r\equiv(dB/dt)/B$). Equation 1 is fairly well reproduced in our simulations, as can be seen from Figure 2$b$, which shows the volume-averaged heating rate of ions $dU_{i}/dt$ (solid-black) and $r\Delta p_{i}$ (solid-green) for run S2m2b0.5.

Since $dU_{i}/dt$ is dominated by the ion pressure anisotropy, the ion heating is ultimately regulated by the pitch-angle scattering provided by either the IC or mirror instabilities. In this section we show that when this scattering is provided mainly by the IC modes, it can also give rise to significant stochastic ion acceleration.

Figure 5$a$ shows the evolution of the ion spectrum for run S2m2b0.5 ($\beta_{i}^{\textrm{init}}=0.5$, $m_{i}/m_{e}=2$), with the color bar indicating the time for each spectrum. This simulation shows the rapid growth of a non-thermal tail that starts once the IC instability grows and saturates ($t\cdot s\sim 1$, as seen from Figure 4$a$). By $t\cdot s=2.5$ the tail can be approximated by a power law $dn_{i}/d\gamma_{i}\propto(\gamma_{i}-1)^{-\alpha_{s}}$ with spectral index $\alpha_{s}\approx 3.4$ plus two bumps ($\gamma_{i}$ is the ion Lorentz factor). The non-thermal tail at $t\cdot s=2.5$ reaches Lorentz factors $\gamma_{i}\sim 10$, and contains $\sim 2-3\%$ of the ions and $\sim 18\%$ of their energy.333After fitting the low energy part of $dn_{i}/d\gamma_{i}$ to a thermal Maxwell-Boltzmann distribution, we define the non-thermal tail through the condition that $dn_{i}/d\gamma_{i}$ is at least a factor $2$ larger than the expectation for the thermal distribution. The solid-black line in Fig. 5$a$ represents the final ($t\cdot s=2.5$) spectra for the analogous run S2m10b0.5 (with $m_{i}/m_{e}=10$ instead of 2). The small difference between the $m_{i}/m_{e}=2$ and 10 cases shows that, as long as the ion parameters $k_{B}T_{i}/m_{i}c^{2}$, $\beta_{i}^{\textrm{init}}$, and $\omega_{c,i}^{\textrm{init}}/s$ are the same, the ion to electron mass ratio does not play a significant role in determining the ion acceleration efficiency. The independence of the acceleration mechanism on $m_{i}/m_{e}$ as well as on $\omega_{c,i}^{\textrm{init}}/s$ will be further tested using 1D simulations in §4.

Figure 5$b$ shows the evolution of the ion spectra for the 2D run S2m2b2 ($m_{i}/m_{e}=2$ and $\beta_{i}^{\textrm{init}}=2$). In this case the growth of the non-thermal energy tail is also present but with a slower growth throughout the whole simulation. By $t\cdot s=2.5$ the tail can be approximated by a much less pronounced power law with spectral index $\alpha_{s}\approx 5.9$. The solid-black line represents the final spectrum for the analogous run S2m10b2, with $m_{i}/m_{e}=10$. The small difference between the $m_{i}/m_{e}=2$ and 10 cases suggests that, as in the IC-dominated case, the ion to electron mass ratio is fairly unimportant in determining the ion acceleration efficiency.

These results strongly suggest that the presence of IC modes is key for the acceleration of ions. In §4 we show that this is indeed the case making use of 1D simulations in which the mirror modes are artificially suppressed. In §4 we will also make use of the low computational cost of 1D simulations to test the effect of using values of $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$ much larger than the ones used in the 2D runs. As we will see, we will find no significant dependence of the ion acceleration on these parameters.

4. 1D shearing simulations

Since the IC modes propagate mainly parallel to the background magnetic field $\langle\textbf{{B}}\rangle$, in this section we study the ion acceleration due to the IC instability by only capturing modes with wave vector k parallel to $\langle\textbf{{B}}\rangle$. We do so by using the 1D version of our shear coordinates simulations (Riquelme et al., 2012). Analogously to what happens in 2D and 3D, in our 1D shearing runs the spatial domain of the simulation evolves with time, following the shearing flow of the plasma. This implies that the domain rotates and stretches so that the resolved k’s are always parallel to $\langle\textbf{{B}}\rangle$. A detailed description of our 1D setup is in Appendix A.

In the next section, we show the suitability of the 1D setup to study problems dominated by modes with $\textbf{{k}}\parallel\textbf{{B}}$ by comparing 1D runs with 2D simulations that are dominated by the IC instability. Then, we use the 1D runs to: $i)$ provide further evidence that the IC modes are the essential ingredient for the ion acceleration, and $ii)$ explore the dependence of the ion acceleration on $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$, which will make use of the low computational cost of the 1D runs.

4.1. 1D vs. 2D comparison

We use runs S1m2b0.5 (1D) and S2m2b0.5 (2D) to compare the 1D and 2D results. First, we check whether our 1D and 2D runs give similar ion spectra. Fig. 6$a$ shows the evolution of the ion spectrum for run S1m2b0.5 (1D) from $t\cdot s=0$ to $t\cdot s=2.5$. We see that the 1D spectral evolution is very similar to the one of the 2D run S2m2b0.5 shown in Fig 5$a$. Indeed, both spectra can be described as a power law of index $\alpha_{s}\approx 3.4$ plus two bumps, with the intermediate energy bump appearing at $t\cdot s\gtrsim 2$. A more detailed comparison can be made by overplotting the final ($t\cdot s=2.5$) spectrum for run S2m2b0.5 in Figure Fig. 6$a$ (black line). We see that the two final spectra are very similar, with the main difference being a $\sim 2$ times larger maximum energy in the 1D run. This small difference is to some extent expected, due to the presence of mirror modes in the 2D case. Indeed, in 2D we have some contributions of mirror modes, which however are not conducive to ion acceleration. This explains why the 2D setup leads to somewhat lower energy gains than in 1D.

In terms of the ion pressure anisotropy, Fig. 2$a$ shows $p_{i,\parallel}$ (red) and $p_{i,\perp}$ (blue) for runs S1m2b0.5 (1D; dotted) and S2m2b0.5 (2D; solid). Both for $p_{i,\parallel}$ and $p_{i,\perp}$, the 1D and 2D simulations give essentially the same results. In dashed lines we show the corresponding double-adiabatic behavior (Chew et al., 1956), which is followed quite well by the two simulations until $t\cdot s\approx 1$. With respect to energy conservation, Figure 2$b$ shows the volume average $dU_{i}/dt$ (black) and $r\Delta p_{i}$ (green) for the same runs S1m2b0.5 (1D; dashed) and S2m2b0.5 (2D; solid). $dU_{i}/dt$ behaves very similarly in 1D and 2D, and in the two cases it corresponds quite well to the heating due to anisotropic viscosity.

In order to compare the behavior of the magnetic fluctuations $\delta\textbf{{B}}$, we use $\delta B_{z}$ since in the 2D runs this component is mainly produced by the IC modes, and it is essentially not affected by the (subdominant but still present) mirror modes. Figure 7$a$ shows a snapshot of $\delta B_{z}$ at $t\cdot s=2$ for the 1D run S1m2b0.5. We see that both in terms of the dominant wavelength ($\approx 30R_{L,i}^{\textrm{init}}$) and of its amplitude, $\delta B_{z}$ behaves fairly similarly to the 2D case, shown in Fig. 3$c$. Figures 7$b$, 7$c$ and 7$d$ show power spectra of $\delta B_{z}$ at different simulation times for a more detailed comparison. In order to reduce the effects of time variability, we take averages during $1.1<t\cdot s<1.2$, $1.8<t\cdot s<2$, and $2.2<t\cdot s<2.4$, respectively. We see that the 1D and 2D spectra look quite similar. Their main difference consists of a small (by a factor $\sim 1.5$) shift in the peak of the 1D spectra towards longer wavelengths, and a factor $\sim 2$ increase in the peak amplitude.

These differences in the $\delta B_{z}$ spectra can be explained to a large extent by the small differences in the ion energy spectra. Indeed, in §5.2 we show that the wave number $k$ at the peak of the IC wave spectrum is determined by the resonance condition with the highest energy ions, with $k\propto 1/\gamma_{i}$. Thus, since the 1D runs produce maximum ion energies $\sim 2$ times larger than in the 2D case, the wave number at the peak should be reduced by a similar factor.

Something similar occurs with the difference in amplitude of $\delta B_{z}$. For relativistic ions interacting resonantly with parallelly propagating waves, the effective scattering frequency $\nu_{eff,i}$ should scale as (Kulsrud & Pearce, 1969):

[TABLE]

Thus, since the ions in the 1D and 2D runs are scattered at roughly the same rate (given their similar evolution of $p_{\parallel,i}$ and $p_{\perp,i}$, as shown in Fig. 2$a$), Eq. 2 implies that the peak value of $dB_{z}^{2}{k}/dk$ in the 1D case should roughly be $\sim 2$ times larger than in the 2D case, which is seen in panels $b$, $c$ and $d$ of Fig. 7.

Thus, besides a factor $\sim 2$ difference in the highest energy of ions (which is likely due to the weak presence of mirror modes in the 2D runs), the 1D runs reproduce reasonably well the 2D results, and provide a valuable tool to study the effect of IC waves on ion acceleration. In the next two sections we use 1D runs to provide further evidence that the IC modes are indeed the essential ingredient for the ion acceleration, and to explore the dependence of the acceleration on $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$.

4.2. The role of IC and mirror modes

Simulations in 1D can be used to further clarify the role of IC and mirror modes in the acceleration of ions. We do this by comparing the ion spectra from a 2D simulation where the mirror modes dominate ($\beta_{i}^{\textrm{init}}=2$) with an analogous 1D run where these modes are artificially suppressed. This is done in Fig. 6$b$, which shows in black the final spectrum of the 2D run S2m2b2, where the mirror modes dominate. In addition, Fig. 6$b$ shows the ion spectra at different times for the 1D run S1m2b2, where the plasma is under the same conditions as in run S2m2b2. It can be seen that ions tend to be significantly more accelerated in the 1D case, in which the IC modes dominate. This shows that the main ingredient for the acceleration of ions is indeed the scattering by the IC modes, with the mirror instability suppressing the acceleration. This suppression is in line with the lack of electric field associated with the mirror modes, as shown in Fig. 3. Thus, when the mirror modes dominate, the scattering of ions tends to be elastic, and no acceleration effect should be present.

4.3. $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}$ dependence

We use 1D simulations to explore the dependence of the ion acceleration by the IC instability on $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$, focusing on the case with $\beta_{i}^{\textrm{init}}=0.5$. In terms of the evolutions of $p_{\perp,i}$ and $p_{\parallel,i}$, Figure 8$a$ shows the cases of runs S1m8b0.5 ($\omega_{c,i}^{\textrm{init}}/s=800$, $m_{i}/m_{e}=8$; dotted line), S1m128b0.5 ($\omega_{c,i}^{\textrm{init}}/s=800$, $m_{i}/m_{e}=128$; solid line) and S1m8b0.5d ($\omega_{c,i}^{\textrm{init}}/s=3200$, $m_{i}/m_{e}=8$; dashed line). No significant difference can be seen between the different mass ratios and magnetizations. The same thing happens when we look at the ion energy gain. Figure 8$b$ shows $dU_{i}/dt$ (black) and $r\Delta p_{i}$ (green) for the same runs. We see that the ion energy gain is fairly independent of $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$, and in all cases it reasonably well agrees with the heating prediction through anisotropic viscosity.

Figure 9$a$ shows the final ion spectra ($t\cdot s=2.5$) for simulations with $m_{i}/m_{e}=2$, 8, 32, and 128 (runs S1m2b0.5, S1m8b0.5, S1m32b0.5, and S1m128b0.5 in Table 1). In all simulations the ions share the same parameters: $\omega_{c,i}^{\textrm{init}}/s=800$, $\beta_{i}^{\textrm{init}}=0.5$, and $k_{B}T_{i}=k_{B}T_{e}=0.05m_{i}c^{2}$, so the only difference is the value of $m_{i}/m_{e}$. The non-thermal ion tail in all cases can be fairly well described as a power-law of spectral index $\alpha_{s}\approx 3.4$ plus two bumps occurring at roughly the same energies. This result shows that, as long as the electrons are somewhat less massive than the ions, their effect on the ion acceleration by IC modes becomes negligible. This is expected given the resonant nature of the interaction between the unstable IC modes and the ions, which requires the modes to be left-handed, circularly polarized (e.g., Gary, 1992). This polarization requirement naturally makes it significantly more difficult for the electrons to interact resonantly with the IC modes, even for a mass ratio as small as $m_{i}/m_{e}=2$.

Similarly, Figure 9$b$ shows the final spectra for simulations with $m_{i}/m_{e}=8$, $\beta_{i}^{\textrm{init}}=0.5$, and $k_{B}T_{i}=k_{B}T_{e}=0.05m_{i}c^{2}$, but with $\omega_{c,i}^{\textrm{init}}/s=400$, 800, 1600 and 3200 (runs S1m8b0.5b, S1m8b0.5, S1m8b0.5c and S1m8b0.5d, respectively). The spectra get slightly harder as $\omega_{c,i}^{\textrm{init}}/s$ increases, with the difference between them being progressively less significant as $\omega_{c,i}^{\textrm{init}}/s$ grows. However, in all cases the tail can be well described as a power-law of spectral index $\alpha_{s}\approx 3.4$ plus two bumps.

The independence of the ion acceleration on $\omega_{c,i}^{\textrm{init}}/s$ can be inferred from the way the effective ion scattering rate $\nu_{eff,i}$ provided by the IC waves is related to $s$. This scattering rate can be estimated from the evolution of $p_{\parallel,i}$ in a homogeneous plasma with no heat flux, assuming that $B$ evolves on time and length scales much larger than $\omega_{c,i}^{-1}$ and $R_{L,i}$ (which is the case in our runs). This evolution is given by Eq. 1 of Sharma et al. (2007):

[TABLE]

where v is the plasma bulk velocity and $\hat{\textbf{{b}}}\equiv\textbf{{B}}/B$. In the case of the shearing plasma ($\textbf{{v}}=-sx\hat{y}$), $\nabla\cdot(p_{\parallel,i}\textbf{{v}})=0$ and $\hat{\textbf{{b}}}\hat{\textbf{{b}}}:\nabla\textbf{{v}}=\hat{b}_{x}\hat{b}_{y}s$. Fig. 8$a$ shows that, after the saturation of the IC modes, $p_{\parallel,i}$ changes at a rate much smaller than $s$ so we can approximate $\partial p_{\parallel,i}/\partial t\approx 0$. Additionally, $\hat{b}_{x}\hat{b}_{y}$ ranges between 0.5 at $t\cdot s=1$ and 0.34 at $t\cdot s=2.5$. Thus, we simply assume $\hat{b}_{x}\hat{b}_{y}\approx 1/2$ and Eq. 3 becomes:

[TABLE]

By comparing the evolutions of $p_{\perp,i}$ and $p_{\parallel,i}$ for runs S1m8b0.8 ($\omega_{c,i}^{\textrm{init}}/s=800$; dot-dashed line) and S1m8b0.8d ($\omega_{c,i}^{\textrm{init}}/s=3200$; dashed line), Figure 8$a$ shows that the factor $p_{\parallel,i}/\Delta p_{i}$ is fairly independent of $\omega_{c,i}^{\textrm{init}}/s$. Thus, Eq. 4 implies that $\nu_{eff,i}\propto s$.

This proportionality between $\nu_{eff,i}$ and $s$ means that the average number of scatterings experienced by the ions after $t=2.5s^{-1}$ (at the end of the simulations) should be about the same in all runs. This property, if the IC modes properties are the same in all the simulations (as it occurs with the runs shown in Fig. 9$b$), should make the accelerating effect of the IC modes independent of $\omega_{c,i}^{\textrm{init}}/s$.

Our 1D simulations, therefore, show that the ion acceleration by IC modes is fairly independent of $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$. In the next section, we use 1D and 2D simulations to describe in further detail the way that ions of different energy get accelerated, emphasizing the role played by their resonant interaction with the IC waves.

5. The acceleration mechanism

The two possible sources of energy for the ions in our simulations are: the energy gain due to anisotropic viscosity (Eq. 1), and the energy gain due to the work done by the electric field associated with the IC waves (considering that the electric field associated to the mirror modes is negligible). In §5.1 we identify the contributions of each of these energy sources to producing the non-thermal ion spectra.

5.1. IC work vs. anisotropic viscosity

Figure 10 shows the different contributions to the energy gain of three ion populations from the 2D run S2m2b0.5, separated according to their final energy at $t\cdot s=2.5$. These populations are:

The “thermal ions”, corresponding to ions in the bulk of the ion distribution, and with their energy gain plotted in Fig. 10$a$. These ions are chosen so that at $t\cdot s=2.5$ their Lorentz factors satisfy $0.09<\gamma_{i}-1<0.1$ (marked by the vertical red line in Fig. 5$a$). 2. 2.

The ions in the “intermediate energy” bump of the tail (shown in Fig. 10$b$). These ions have Lorentz factors in the range $1<\gamma_{i}-1<1.03$ at $t\cdot s=2.5$ (marked by the vertical dark green line in Fig. 5$a$). 3. 3.

The ions in the “high energy” bump of the tail (shown in Fig. 10$c$). This population corresponds to the highest energy ions, defined by $\gamma_{i}-1>4.8$ at $t\cdot s=2.5$ (marked by the grey region in Fig. 5$a$).

For each of these populations we plot the following contributions to their energy gain:

The work done by the electric field of the IC waves, $W_{IC}$, which is shown by the blue lines of Fig. 10. 2. 2.

The energy gain by anisotropic viscosity, $AV$, shown by the green line in Fig. 10.444For each ion population, this energy gain is calculated as the integral in time of the rate of energy gain due to viscosity: $\int dt\,r\,\Delta p_{i}$ (see Eq. 1), where $\Delta p_{i}$ is calculated using only the ions of each population.

The blue line in Figure 10$a$ shows that the energy given by the electric field of the IC waves to the thermal ions, $W_{IC}$, is negative. This implies that the scattering process on average substracts energy from the thermal ions and transfers it to the waves. The total gain in energy of the thermal ions is still positive, and dominated by viscous heating. On the other hand, Figure 10$b$ shows that the work done by the IC waves on the ions of the intermediate energy bump is positive and larger than the heating by anisotropic viscosity, $AV$, which means that these ions are mainly energized by the scattering caused by the IC waves. This energization occurs mainly after $t\cdot s\approx 2$, which is consistent with the late time appearance of the intermediate energy bump, as shown by the time evolution of the ion spectrum depicted in Figure 5$a$. Finally, Figure 10$c$ shows that $W_{IC}$ is about three times larger than $AV$, implying that the IC acceleration for the highest energy ions is even larger than for the “intermediate energy” ions.

In the three Figures 10$a$, 10$b$ and 10$c$ we also plot the sum of $AV$ and $W_{IC}$ (red line) and the average change in energy of the three ion populations, $\langle\Delta E_{i}\rangle$ (black line). We see that these two quantities are essentially the same for the three populations, implying that the energy gain due to anisotropic viscosity and the electric field of the IC modes accounts quite well for the total ion energy evolution in the three populations.

5.2. Resonance with IC waves

In §5.1 we show that a non-thermal ion tail is produced by the scattering of ions off IC waves, which, in turn, obtain their energy from the pressure anisotropy of the thermal ions. This implies that this acceleration mechanism requires the resonance condition between ions and IC waves to be satisfied by both the thermal and non-thermal ions. The resonance condition is:

[TABLE]

where $\omega$ and $k$ are the real part of the frequency and the wave number of the modes, $v_{||}$ is the ion velocity parallel to B, and $\omega_{c,i}$ ($\equiv eB/m_{i}c$) is the non-relativistic cyclotron frequency of the ions.

In order to check that this resonance condition is satisfied by ions of all energies, in Figs. 11$a$ and 11$b$ we measure the ratio $\omega/k$ for the IC waves at two different time intervals for the 1D run S1m2b0.5 ($\beta_{i}^{\textrm{init}}=0.5$, $m_{i}/m_{e}=2$, $\omega_{c,i}^{\textrm{init}}=800$). We do this by plotting $|\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)|^{2}$, where $\delta B_{z}(\omega,k)$ and $\delta B_{\perp,xy}$ correspond to two mutually perpendicular components of $\delta\textbf{{B}}$, that are also perpendicular to $\langle\textbf{{B}}\rangle$ (see caption of Fig. 4), the tilde ($\,\tilde{}\,$) denotes the time and space Fourier transform of a quantity and $i=\sqrt{-1}$. Thus, the combination $\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)$ allows to separate the contributions to $\delta\textbf{{B}}$ provided by IC waves (left-handed, circularly polarized) and whistler waves (right-handed, circularly polarized), with the latter being expected to be destabilized by the pressure anisotropy of electrons (Gary & Wang, 1996). In the case of run S1m2b0.5, the IC waves contribute to $|\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)|$ only for $\omega>0$, while the whistler wave contribution appears for $\omega<0$. This way, calculating $|\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)|^{2}$ allows to separate the IC and whistler contributions to $\delta\textbf{{B}}$, and to estimate $\omega(k)$ for these two modes.

Figure 11$a$ corresponds to the time interval $1.1<t\cdot s<1.2$ of run S1m2b0.5. The IC modes have a phase velocity of $\omega/k\approx 0.5c$. Additionally, we measure the rms ion velocity parallel to B at $1.1<t\cdot s<1.2$, which is $v_{\parallel}^{\textrm{rms}}\approx 0.16c$,555This can be estimated by the factor $\sim 2$ decrease in $p_{i,\parallel}$ seen in Fig. 2$a$, and considering that initially $k_{B}T_{i}/m_{i}c^{2}=0.05$. implying that to a good approximation we can neglect the $v_{\parallel}$ term on the left hand side of Eq. 5. The resonance condition at $1.1<t\cdot s<1.2$ can thus be written as:

[TABLE]

At $1.1<t\cdot s<1.2$, $\gamma_{i}$ is in the range $1<\gamma_{i}\lesssim 2$ (see Fig. 6$a$), implying that most of the power of the IC modes should be in the range $1\lesssim kc/\omega_{c,i}\lesssim 2$, which coincides well with the range of $k$ in which most of the IC power is observed in Fig. 11$a$.

Analogously, Fig. 11$b$ shows $|\delta B_{z}(\omega,k)+i\delta B_{\perp,xy}(\omega,k)|^{2}$ for the same simulation but in the time range $2.2<t\cdot s<2.3$. In this case, Fig. 11$b$ shows that $\omega/k\approx 0.7c$. Thus, making a similar analysis as in the case $1.1<t\cdot s<1.2$, we obtain:666Here we also assume $\omega/k$ ($\approx 0.7c$) $\gg v_{\parallel}$, which allows neglecting $v_{\parallel}$ in Eq. 5. This is a reasonable approximation considering that at $2.2<t\cdot s<2.3$, $v_{\parallel}$ (measured directly from the simulation) is always smaller than $\sim 0.4c$, even considering the highest energy particles in the non-thermal tail.

[TABLE]

Since in this time interval $1<\gamma_{i}\lesssim 8$, we obtain that most of the power of the IC modes should be in the range $0.2\lesssim kc/\omega_{c,i}\lesssim 1.5$. This interval coincides reasonably well with the range of $k$ where most power is concentrated in Fig. 11$b$. Notice that this power appears to be enhanced in two intervals of $k$. The high-$k$ interval corresponds to $0.7\lesssim kc/\omega_{c,i}\lesssim 1.5$, which, according to Eq. 7, implies resonance with ions with $1<\gamma_{i}\lesssim 2$. Remarkably, this is a range of $\gamma_{i}$ with abundant IC scattering at $2.2\lesssim t\cdot s\lesssim 2.3$, as shown by the rapid formation of the “intermediate energy” bump, which mainly occurs at $2\lesssim t\cdot s\lesssim 2.5$. The low-$k$ enhancement occurs for $0.2\lesssim kc/\omega_{c,i}\lesssim 0.5$, which, according to Eq. 7, corresponds to $3\lesssim\gamma_{i}\lesssim 8$. This $\gamma_{i}$ interval coincides well with the “high energy” bump shown in Fig. 6$a$ at $t\cdot s\sim 2.2-2.3$, and is also consistent with the rapid increase in energy of this bump.

Thus, we have shown that the range of $k$ in which the amplitude of the IC modes is significant is consistent with the resonance condition occurring for both the thermal and non-thermal ions in the tail. Figs 11$c$ and 11$d$ show the same quantities as Figs 11$a$ and 11$b$ but for simulation S1m8b0.5, where the ions are under the same conditions as in run S1m2b0.5, but with $m_{i}/m_{e}=8$ instead of $m_{i}/m_{e}=2$. The quantity $|\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)|^{2}$ essentially preserves the same properties for the IC modes ($\omega>0$). The fact that $\omega(k)$ of the IC modes is fairly independent of $m_{i}/m_{e}$ is consistent with the near independence of the ion acceleration process on $m_{i}/m_{e}$. The whistler modes, on the other hand, do change significantly their properties as $m_{i}/m_{e}$ is increased, for instance by getting weaker and increasing their frequency. This is expected since quantities like $\omega_{c,e}/\omega_{p,e}$, $\omega_{c,e}/s$ and $k_{B}T_{e}/m_{e}c^{2}$ do change when varying $m_{i}/m_{e}$.

Finally, the four panels in Fig. 11 also include IC and whistler theoretical dispersion relation calculations, $\omega_{\textrm{theo}}(k)$ (dashed lines), obtained with the linear Vlasov solver NHDS (Verscharen & Chandran, 2018). These calculations assume the ion and electron conditions obtained in the simulations in terms of their temperatures, pressure anisotropies, mass ratio and Alfvén velocity. However, they do not consider departures from Maxwell-Boltzmann distributions or relativistic effects. In the four panels $\omega_{\textrm{theo}}(k)$ reproduces the behavior of $|\delta\tilde{B}_{z}(\omega,k)+i\delta\tilde{B}_{\perp,xy}(\omega,k)|^{2}$ for the case of the IC modes well, which shows that the phase (and group) velocity of the IC waves obtained from the simulations are not significantly affected by the non-thermal or relativistic effects in the ion velocity distribution. Most discrepancies occur for whistler waves when $m_{i}/m_{e}=8$. In this case, the discrepancy is most likely due to the electrons being significantly relativistic (Lorentz factor $\gamma_{e}\gg 1$), which is a regime strictly not captured by NHDS.

6. Comparing shear vs. compression

The ion acceleration presented in this paper occurs during the non-linear, saturated stage of the IC instability. This stage is reached through the continuous amplification of a background magnetic field $\langle\textbf{{B}}\rangle$, which in previous sections has been driven by an imposed shear plasma motion. In this section we show that, as long as the ion conditions are similar, the specific process that amplifies $\langle\textbf{{B}}\rangle$ does not play an important role in the acceleration. We do this by running 1D simulations in which $\langle\textbf{{B}}\rangle$ grows due to plasma compression. The simulation setup is the same as in Sironi & Narayan (2015). The plasma is compressed along the $y$ and $z$ axes, with $\langle\textbf{{B}}\rangle$ pointing along $x$. In this setup, $\langle\textbf{{B}}\rangle$ is amplified due to magnetic flux conservation (see Fig. 1 of Sironi & Narayan, 2015). In this situation, the background field grows as $|\langle\textbf{{B}}\rangle|=B_{0}(1+qt)^{2}$, where $q$ is a constant that roughly corresponds to the compression rate of the plasma. These 1D simulations only capture modes with $\textbf{{k}}\parallel\langle\textbf{{B}}\rangle$, which means that the mirror modes are artificially suppressed.

The compression runs use the same initial plasma parameters as the shearing runs that show significant ion acceleration: $\beta_{i}^{\textrm{init}}=0.5$, $k_{B}T_{i}/m_{i}c^{2}=0.05$, and $T_{i}=T_{e}$, which allows a direct comparison between the two setups. Fig. 12$a$ shows the evolution of $T_{\perp,i}$ (solid-blue) and $T_{\parallel,i}$ (solid-red) for run C1m8b0.5a, which uses $m_{i}/m_{e}=8$ and $\omega_{c,i}^{\textrm{init}}/q=1600$ ($k_{B}T_{\perp,i}=p_{\perp,i}/n_{i}$ and $k_{B}T_{\parallel,i}=p_{\parallel,i}/n_{i}$, where $n_{i}$ is the ion density).777We show $T_{\perp,i}$ and $T_{\parallel,i}$ instead of $p_{\perp,i}$ and $p_{\parallel,i}$ in order to disentangle variations in $n_{i}$ through the compression from the variation through double-adiabatic and instability-induced effects. Both temperatures initially follow the double adiabatic evolution reasonably well, which is marked with the dashed-blue and dashed-red lines, respectively (Chew et al., 1956).888The small initial discrepancies ($t\cdot s\lesssim 0.6$) between the double-adiabatic predictions and $T_{\parallel,i}$ and $T_{\perp,i}$ are because these predictions assume that the ions are non-relativistic, which is mildly broken for our $k_{B}T_{i}/m_{i}c^{2}\sim 0.1$ ions. The adiabatic evolution breaks at $t\cdot q\approx 0.6$, after the exponential growth of the IC modes begins. This exponential growth and subsequent saturation can be seen from Fig. 12$b$, which shows the $\delta B_{y}^{2}$ and $\delta B_{z}^{2}$ components of run C1m8b0.5a, which evolve quite similarly to the $\delta B_{\perp,xy}$ and $\delta B_{z}^{2}$ fluctuations of analogous shearing runs (see, e.g., the cases of runs S2m2b0.5 and S2m10b0.5 in Fig. 4$a$).

The evolution of the ion spectrum of run C1m8b0.5a is shown in Fig. 12$c$. We see a remarkable similarity with the spectral evolution of the shearing runs with $\beta_{i}^{\textrm{init}}=0.5$, shown in Figs. 5$a$ and 6$a$. Indeed, also in the compressing case, there is a growth of a non-thermal tail that, by $t\cdot q=2.5$, can be described as a power law of index $\alpha\approx 3.5$ plus two bumps.

We also test the dependence of this acceleration on both $m_{i}/m_{e}$ and ion magnetization $\omega_{c,i}^{\textrm{init}}/q$. Fig. 12$d$ shows a comparison between cases with $m_{i}/m_{e}=8$ and 16 (black and red, respectively), and with $\omega_{c,i}/q=1600$ and 3200 (solid and dashed, respectively), which share the same initial parameters: $\beta_{i}^{\textrm{init}}=0.5$, $k_{B}T_{i}/m_{i}c^{2}=0.05$ and $T_{i}=T_{e}$. We see essentially no difference between the runs with different values of $m_{i}/m_{e}$, and only a slight hardening of the tail as $\omega_{c,i}/q$ increases, in agreement with the shearing results of §4.3.

These results essentially reproduce our shearing runs, both in terms of the significance of the ion acceleration in the regime: $\beta_{i}^{\textrm{init}}=0.5$, $k_{B}T_{i}/m_{i}c^{2}=0.05$ and $T_{i}=T_{e}$, and in terms of the almost independence on $m_{i}/m_{e}$ and $\omega_{c,i}/q$.

As in the shearing case, the independence on $m_{i}/m_{e}$ can be understood as due to the resonant nature of the ion-IC interaction. The independence on $\omega_{c,i}^{\textrm{init}}/q$ implies, similarly to the shearing case, that the effective ion scattering rate $\nu_{eff,i}$ has to be $\propto q$. This can indeed be inferred by applying Eq. 3 to the compressing runs. In this case, $\hat{\textbf{{b}}}\hat{\textbf{{b}}}:\nabla\textbf{{v}}=0$, thus (using the continuity relation $\partial n_{i}/\partial t=-n_{i}\nabla\cdot\textbf{{v}}$) $n_{i}\partial T_{\parallel,i}/\partial t=(2/3)\nu_{eff,i}\Delta p_{i}$. However, Fig. 12$a$ shows that, after the saturation of the IC waves, $\partial T_{\parallel,i}/\partial t\sim(q/2)T_{\parallel,i}$, which suggests:

[TABLE]

The ratio $p_{\parallel,i}/\Delta p_{i}$ evolves fairly independent of $\omega_{c,i}^{\textrm{init}}/q$, as can be seen from Fig. 12$a$, which shows $T_{\perp,i}$ and $T_{\parallel,i}$ for runs C1m8b0.5a ($\omega_{c,i}^{\textrm{init}}/q=1600$; solid line) and C1m8b0.5b ($\omega_{c,i}^{\textrm{init}}/q=3200$; dot-dashed line). Therefore, Eq. 8 implies that $\nu_{eff,i}\propto q$.

The similarity between Eqs. 4 and 8, is consistent with the ion spectra being similar in the two setups, when comparing spectra at equal values of $ts$ and $tq$. Indeed, considering that in the shearing and compressing cases $p_{\parallel,i}/\Delta p_{i}\sim 1/2$ during most of the saturated IC regime (see Figs. 2$a$ and 12$a$), at equal values of $ts$ and $tq$, the ions must have experienced a similar number of effective scatterings. Thus, if the properties of the IC modes are comparable, their acceleration effects by the end of the shearing and compressing simulations should also be comparable.999Notice, however, that, although the initial ion conditions in the shearing and compressing runs are the same, the final conditions are somewhat different. For instance, in the shearing runs with $\beta_{i}^{\textrm{init}}=0.5$ the final value of the parallel ion beta is $\beta_{\parallel,i}^{\textrm{final}}\sim 0.025$ (considering the evolution of $p_{\parallel,i}$ seen in Fig. 2$a$ and the expected evolution of $B$), while in the analogous compressing runs $\beta_{\parallel,i}^{\textrm{final}}\sim 0.1$ (considering the evolution of $T_{\parallel,i}$ seen in Fig. 12$a$ and the expected evolutions of $n_{i}$ and $B$). This, plus the different factors on the right hand sides of Eqs. 4 and 8, imply that the final ion spectra in these two setups should be similar but not necessarily the same.

7. Summary and Conclusions

Our 1D and 2D particle-in-cell (PIC) plasma simulations show that ions can be stochastically accelerated by the inelastic scattering provided by the ion-cyclotron (IC) instability. This acceleration occurs in the non-linear, saturated state of the instability, which is reached due to a permanent amplification of the background magnetic field $\langle\textbf{{B}}\rangle$.

In the regime in which initially $k_{B}T_{i}=k_{B}T_{e}=0.05m_{i}c^{2}$, we show that the IC ion acceleration is significant if $\beta_{i}^{\textrm{init}}\lesssim 1$. This is demonstrated by comparing 2D simulations with $\beta_{i}^{\textrm{init}}=0.5$ and 2. When $\beta_{i}^{\textrm{init}}=0.5$, the ion scattering is dominated by the IC instability, which produces a non-thermal tail in the ion energy spectrum. After B is amplified by a factor $\sim 2.7$, the tail can be approximately described as a power-law of index $\sim 3.4$ plus two non-thermal bumps. The maximum ion Lorentz factor at that time is $\gamma_{i}\sim 10$, but it continues to grow at the end of the simulation. Also, the tail accounts for $2-3\%$ of the ions and $\sim 18\%$ of their kinetic energy. On the other hand, when $\beta_{i}^{\textrm{init}}=2$, the ion scattering is dominated by the mirror instability (the IC modes are subdominant) and the acceleration is significantly suppressed.

In the IC dominated regime, as the ion scattering increases the energy of the ions of the tail, it reduces the energy of the ions in the thermal part of the spectrum (see Fig. 10). This is consistent with the IC modes being driven unstable mainly by the pressure anisotropy of the thermal ions. This way, the role of the IC modes is to absorb part of the energy of the thermal ions and give it to the non-thermal ions. This process is very similar to the stochastic acceleration of electrons by the whistler instability found by Riquelme et al. (2017). The efficiency of the ion acceleration, therefore, relies on the IC modes being able to provide resonant scattering to both thermal and non-thermal ions. We analyzed the consistency of this scenario by calculating the $k$ numbers and phase velocities of the dominant IC modes, showing that they can resonate with ions of all the energies (see discussion in §5.2).

Given that our simulations can not use realistic values of $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$, one important aspect of our study is to ensure that these parameters do not affect the acceleration. Thus, first we ensure that the dominance of the IC modes for $\beta_{i}^{\textrm{init}}=0.5$ does not depend on $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$. This was done comparing 2D simulations with $m_{i}/m_{e}=2$ and 10 (see §3.1), and also using theoretical, linear dispersion relation calculations to determine the pressure anisotropy needed for the growth of the IC and mirror instabilities in astrophysically realistic conditions (see Appendix §B). Both analyses show that, in realistic astrophysical plasmas, the IC instability dominates in the regime $\beta_{i}^{\textrm{init}}\lesssim 1$, at least for the case of $T_{e}=T_{i}$ explored here.

Then, using that the dominant IC wave vectors k satisfy $\textbf{{k}}\parallel\langle\textbf{{B}}\rangle$, we use computationally cheaper 1D shear simulations to test the ion acceleration using a significantly larger range of values for $m_{i}/m_{e}$ and $\omega_{c,i}^{\textrm{init}}/s$. While $m_{i}/m_{e}$ almost does not affect the acceleration, increasing $\omega_{c,i}^{\textrm{init}}/s$ only produces a slight hardening of the non-thermal tail (see §4.3). This almost complete independence of the acceleration on $\omega_{c,i}^{\textrm{init}}/s$ is consistent with the effective ion scattering rate $\nu_{eff,i}$ being proportional to $s$. This condition is indeed needed in order to have the continuous driving of the ion pressure anisotropy being nearly compensated by the pitch-angle scattering.

In order to assess the importance of the specific large scale mechanism that amplifies the background magnetic field, we also ran compressing box PIC simulations like in Sironi & Narayan (2015). We find essentially no difference in the ion acceleration efficiency between the shearing and compressing cases (see §6).

Our work is valid in a sub-relativistic regime in which initially $k_{B}T_{i}=k_{B}T_{e}=0.05m_{i}c^{2}$. This regime can be relevant in the inner region of low-luminosity accretion disks around black holes (where the collisionless plasma condition is expected). In these systems, the condition $\beta_{i}^{\textrm{init}}\lesssim 1$ required for the acceleration is most likely satisfied in the coronal region of the disks (e.g., Chael et al., 2018).

Nevertheless, assessing the importance of the presented acceleration mechanism in these and other astrophysical systems requires a more complete understanding of its dependence on plasma parameters, as well as clarifying the importance of possible long term evolution effects. Indeed, in this work we focus on a single value of $k_{B}T_{i}$ and use $T_{e}/T_{i}=1$. However, varying these parameters may affect the IC physics and, therefore, the efficiency of the ion acceleration. For instance, having $T_{e}/T_{i}\ll 1$ may increase significantly the values of $\beta_{i}^{\textrm{init}}$ for which the IC instability dominates (Sironi & Narayan, 2015; Sironi, 2015). Since the condition $T_{e}/T_{i}\ll 1$ is most likely satisfied in low luminosity disks (Narayan & Yi, 1995; Yuan et al., 2003), this could increase the importance of the presented ion acceleration in these systems.

In terms of the long term evolution of the acceleration process, in the turbulent environment of accretion disks we expect many successive $\delta B/B\sim 1$ amplifications and decreases of the field (see also Verscharen at al., 2016). So a more realistic picture of this process should consider the acceleration presented in this work occurring many times as the plasma is gradually accreted. We will study these aspects of the acceleration process in future investigations.

This research was supported by the supercomputing infrastructure of the NLHPC (ECM-02) at the Center for Mathematical Modeling of University of Chile, by the Habanero cluster at Columbia University, and by the XSEDE computing system (allocations TG-AST140039, TG-AST140083 and TG-PHY160040). F.L. acknowledges support from the NSF grant NST AST-1616037. L.S. thanks support from DoE DE-SC0016542, NASA ATP NNX-17AG21G, NSF ACI1657507, and NSF AST-1716567. D.V. was supported by the STFC Ernest Rutherford Fellowship ST/P/003826/1.

Appendix A A. 1D shear setup

The goal of our 1D shear runs is to simulate a shearing plasma assuming that its properties depend only on the direction parallel to the mean background magnetic field $\langle\textbf{{B}}\rangle$. This is equivalent to assuming that the wave vectors k captured in the simulations satisfy $\textbf{{k}}\parallel\langle\textbf{{B}}\rangle$. Since in a shearing plasma the direction of $\langle\textbf{{B}}\rangle$ evolves with time, the orientation of the wave vectors that can be consistently resolved needs to evolve accordingly.

Our 1D setup is built upon the “shearing coordinates” setup presented by Riquelme et al. (2012). This setup was designed so that the simulation domain follows the shearing flow of the plasma, which is given by the shear velocity v ($=v\hat{y}$). Figures 1$a$ and 1$b$ illustrate the way the shape of the domain evolves as seen by an inertial observer at $t=0$ and $t>0$, respectively. Formally, the shearing coordinates $(x^{\prime},y^{\prime},z^{\prime},t^{\prime})$ are defined in terms of the regular, inertial coordinates $(x,y,z,t)$ as

[TABLE]

where $\Gamma=(1-v^{2}/c^{2})^{-1/2}$, $v=-c\textrm{ arctanh}(sx/c)$,101010Notice that, in the limit $|v|\ll c$, the $v=-c\textrm{ arctanh}(sx/c)$ expression is equivalent to $v=-sx$, as we assume in Figures 1$a$ and 1$b$. Since $v=-c\textrm{ arctanh}(sx/c)$ ensures that $|v|<c$, this expression was adopted in Riquelme et al. (2012) in order to deduce the equations that describe the dynamics of the plasma in the shearing coordinate system. However, both in Riquelme et al. (2012) and in this paper we are interested in the plasma dynamics in the regime $|v|\ll c$. $c$ is the speed of light and $s$ is the shear rate of the plasma.

Our 1D domain is defined by $y^{\prime}=0$, and it is shown by the blue lines in Figures 1$a$ and 1$b$. These figures show how the length and orientation of the 1D domain change over time. This in turn changes the orientation of the wave vectors k that can be captured within the domain. In addition, there is an initially homogeneous magnetic field, B, pointing parallel to the 1D domain at $t=0$. Magnetic flux conservation ensures that $\langle\textbf{{B}}\rangle$ will always be parallel to the 1D domain at $t>0$, as also depicted in Figures 1$a$ and 1$b$. Therefore, if the problem of interest is dominated by waves that propagate parallel to $\langle\textbf{{B}}\rangle$, our 1D simulations will be able to capture the essence of the phenomenon.

In order for our 1D simulations to consistently satisfy $\textbf{{k}}\parallel\langle\textbf{{B}}\rangle$, we need to replace $x^{\prime}$ with a new coordinate $x^{\prime}_{1}$, which we define as

[TABLE]

Indeed, if we assume $|v|\ll c$ and $sy^{\prime}\ll c$, the partial derivatives with respect to the spatial coordinates in the inertial frame are

[TABLE]

Thus, if we impose the 1D condition, namely that the fields depend only on $x^{\prime}_{1}$, with $\partial/\partial y^{\prime}=\partial/\partial z^{\prime}=0$, Equations A3 imply that the gradient of any field component in the inertial frame will have coordinates proportional to the vector ($1,-st,0$). In the inertial frame the background magnetic field $\langle\textbf{{B}}\rangle$ has components ($B_{0},0,0$) at $t=0$, then, due to magnetic flux freezing, at $t>0$ these components will be $B_{0}$($1,-st,0$). Thus, after imposing $\partial/\partial y^{\prime}=\partial/\partial z^{\prime}=0$, our 1D simulations are able to capture all the modes with wave vectors k parallel to $\langle\textbf{{B}}\rangle$, despite the fact that $\langle\textbf{{B}}\rangle$ changes orientation over time.

The assumptions $|v|\ll c$ and $sy^{\prime}\ll c$ are equivalent to assuming that the plasma region of interest has a typical size, $L$, that satisfies $L\ll c/s$. For the study of kinetic instabilities, $L$ is typically of the order of the Larmor radius of the particles. Thus, the restriction $L\ll c/s$ applied to species $j$ becomes: $[\gamma_{j}s/\omega_{c,j}][v_{j,\perp}/c]\ll 1$, where $\gamma_{j}$, $\omega_{c,j}$, and $v_{j,\perp}$ are the Lorentz factor, (non-relativistic) cyclotron frequency, and velocity perpendicular to B of particles $j$, respectively. In typical astrophysical environments, $\omega_{c,j}$ is many orders of magnitude larger than $s$, therefore the assumption $L\ll c/s$ is reasonable for the study of kinetic instabilities.

In what follows we will give expressions for the evolution of fields and particle momenta and positions in our 1D shear setup. The condition $L\ll c/s$ is assumed in all of these expressions.

A.1. Evolution of $\textbf{{E}}^{\prime}$ and $\textbf{{B}}^{\prime}$

Considering the evolutions of $\textbf{{E}}^{\prime}$ and $\textbf{{B}}^{\prime}$ in terms of $\textbf{{J}}^{\prime}$ in the shearing coordinates presented in Equations A14 and A26 of Riquelme et al. (2012), after replacing $x^{\prime}$ with $x^{\prime}_{1}$ (Eq. A2) and assuming $\partial/\partial y^{\prime}=\partial/\partial z^{\prime}=0$, the fields dynamics is given by:

[TABLE]

and

[TABLE]

In Equations A4 and A5, $\textbf{{E}}^{\prime}$, $\textbf{{B}}^{\prime}$, and $\textbf{{J}}^{\prime}$, are defined according to the conventional relativistic transformation of the electric field, magnetic field, and current density from the inertial frame to the frame of an observer that moves with velocity v (Eqs. A2 and A28 of Riquelme et al. (2012)). This means that $\textbf{{J}}^{\prime}$ can be calculated directly from the motion of the particles in the shear coordinates, assuming that their momenta $\textbf{{p}}^{\prime}$ can be obtained from their momenta in the inertial frame p, using the conventional relativistic momentum transformation (see Eq. A30 of Riquelme et al. (2012)).

A.2. Evolution of particle momenta

Since our 1D setup uses the same definitions of time ($t^{\prime}$) and of particle momenta ($\textbf{{p}}^{\prime}$) as in the shearing coordinate system of Riquelme et al. (2012), the evolution of $\textbf{{p}}^{\prime}$ will not change and will be given by

[TABLE]

where $q_{c}$ is the particles’ electric charge and $\textbf{{u}}^{\prime}=\textbf{{p}}^{\prime}/m\gamma^{\prime}$, with $m$ and $\gamma^{\prime}$ being the particle mass and Lorentz factor, respectively. Eq. A6 corresponds to Eqs. A31 of Riquelme et al. (2012) in the limit $L\ll c/s$.

A.3. Evolution of particle positions

By taking the derivative of Equation A2 with respect to $t^{\prime}$, and applying the limit $L\ll c/s$ to the expressions for $dx^{\prime}/dt^{\prime}$ and $dy^{\prime}/dt^{\prime}$ given by Eq. A35 of Riquelme et al. (2012), we find

[TABLE]

where $u_{x}^{\prime}$ and $u_{y}^{\prime}$ are the $x$ and $y$ components of $\textbf{{u}}^{\prime}$. Eq. A7 implies that $dx^{\prime}_{1}/dt^{\prime}$ is the scalar product between $\textbf{{u}}^{\prime}$ and $\langle\textbf{{B}}\rangle/|\langle\textbf{{B}}\rangle|$ ($=(\hat{x}-st^{\prime}\hat{y})/(1+s^{2}t^{\prime 2})^{1/2}$), corrected by the “expansion” of the 1D domain (represented by the $y^{\prime}=0$ region depicted in blue in Figures 1$a$ and 1$b$), which provides the extra common factor $1/(1+s^{2}t^{\prime 2})^{1/2}$ on the right hand side of Eq. A7. This implies that $x^{\prime}_{1}$ correctly represents particle displacements along $\langle\textbf{{B}}\rangle$.

A.4. Charge conservation

The evolution of particle positions in our simulations is consistent with the conservation of electric charge density. This can be checked by considering Gauss’ law in the inertial frame: $\nabla\cdot\textbf{{E}}=4\pi\rho_{c}$, where E and $\rho_{c}$ are the electric field and charge density in that frame. Indeed, applying the transformations given by Eqs. A1 to Gauss’ law in the limit $L\ll c/s$, and assuming first the more general problem in 3D described in Riquelme et al. (2012):

[TABLE]

where $\rho_{c}^{\prime}$ is the charge density of the plasma calculated using the shearing coordinates, which, as shown in Eq. A29 of Riquelme et al. (2012), is equal to $\rho_{c}$ in the limit $L\ll c/s$. Thus, according to Eq. A8, in order to satisfy charge conservation in our multidimensional shearing coordinate, particle motions in $y^{\prime}$ lead to an effective current along the $y^{\prime}$ direction of $J^{\prime}_{y}+st^{\prime}J^{\prime}_{x}$. Thus, charge conservation is ensured in our multidimensional shearing coordinate simulations by evolving $x^{\prime},y^{\prime}$ and $z^{\prime}$ according to

[TABLE]

which we do in the 2D simulations presented in this paper and in our previous works (e.g., Riquelme et al., 2017). Eq. A9 is equivalent to Eq. A35 of Riquelme et al. (2012) in the limit $L\ll c/s$.111111Notice that the notation used in Eq. A35 of Riquelme et al. (2012) is different from the one used in this work, and that we interprete their $u_{x}^{\prime},u_{y}^{\prime}$ and $u_{z}^{\prime}$ as our $dx^{\prime}/dt^{\prime}$, $dy^{\prime}/dt^{\prime}$ and $dz^{\prime}/dt^{\prime}$.

Finally, by making an analogous analysis in the case of the 1D setup, we show that charge conservation implies

[TABLE]

Thus, according to Eq. A10, in order for our 1D runs to conserve charge, $x^{\prime}_{1}$ must evolve so that the effective current along $\langle\textbf{{B}}\rangle$ is $(J_{x}^{\prime}-st^{\prime}J_{y}^{\prime})/(1+s^{2}t^{\prime 2})$, which is indeed satisfied by the evolution of $x_{1}^{\prime}$ given by Eq. A7.

Appendix B B. IC vs. mirror under realistic conditions: linear theory analysis

In order to inquire whether the dominance of the IC instability for $\beta_{i}^{\textrm{init}}\lesssim 1$ continues to be valid in realistic, astrophysically relevant regimes ($m_{i}/m_{e}=1836$ and $\omega_{c,i}/s\gg 800$, which we can not study with our 2D simulations), we make use of linear theory. Thus, we calculate the threshold ion pressure anisotropy, $\Delta p_{i}/p_{\parallel,i}\equiv(p_{\perp,i}-p_{\parallel,i})/p_{\parallel,i}$, needed for the growth of IC and mirror modes at a rate $s$, using the linear Vlasov solver NHDS (Verscharen & Chandran, 2018, which assumes bi-Maxwellian, non-relativistic ion velocity distributions).121212Modern linear solvers to the linear Vlasov-Maxwell system of equations can account for relativistic effects and non-Maxwellian background distributions (Verscharen et al., 2018). Studies using this type of solver can evaluate these effects; however, their application is beyond the scope of this paper. This condition is motivated by the assumption that, in order to maintain $\Delta p_{i}/p_{\parallel,i}$ at a nearly stationary level (as we see in Figs. 13$a$ and 13$b$), the modes that provide the pitch-angle scattering must grow at roughly the rate at which the anisotropy is driven, which is $\sim s$. Thus, if non-linear effects did not play any significant role, the dominant instability should be the one with the lowest theoretical anisotropy threshold for a given value of $s$.

However, non-linear effects are expected to be important and to affect the IC and mirror instabilities differently. Indeed, in-situ observations in the solar wind show significant discrepancies between linear theory and the measured ion anisotropy in regions of parameter space in which the IC instability should dominate (Hellinger et al., 2006; Bale et al., 2009). One possible explanation is the departure from bi-Maxwellian ion velocity distributions observed in the solar wind, which may affect the efficiency of the resonant scattering between ions and the IC modes (Isenberg et al., 2013). Another possibility is the inhomogeneity in the magnetic field produced by the (subdominant but still present) mirror modes, which may also affect this resonance (Southwood & Kivelson, 1993). Our approach is thus to estimate these non-linear effects using the simulated cases. These simulations thus provide us with a calibration of the linear theory criterion for determining the dominant instability, which can then be applied to astrophysically realistic regimes.

Figures 13$a$ and 13$b$ show the linear anisotropy threshold given by the mirror (solid red) and IC (solid blue) instabilities with growth rate $s$ in runs S2m2b0.5 and S2m2b2, which have $\beta_{i}^{\textrm{init}}=0.5$ and 2, respectively (both with $m_{i}/m_{e}=2$ and $\omega_{c,i}^{\textrm{init}}/s=800$), and compare them with the anisotropies obtained from the simulations (solid green). First, in the IC dominated regime (run S2m2b0.5), the IC threshold is $\sim 3$ times smaller than the ion anisotropy obtained from the simulation, showing that, in the case of the IC instability, non-linear effects give rise to ion anisotropies significantly larger than what is implied by the linear theory threshold. Also, the IC threshold is at least $\sim 2$ times smaller than the mirror threshold (Fig. 13$a$), while in the mirror dominated case (Fig. 13$b$) the IC threshold is at most $\sim 1.5$ times smaller. These results suggest that, in order for the IC instability to dominate, the IC threshold should be at least $\sim 2$ times smaller than the mirror threshold.

We thus apply this criterion to astrophysically realistic cases, in which $m_{i}/m_{e}=1836$ and the instabilities grow at a rate of $10^{-6}\omega_{c,i}$. Figures 13$a$ and 13$b$ show in dotted lines the corresponding linear anisotropy thresholds for the mirror (dotted red) and IC (dotted blue) instabilities for $\beta_{i}^{\textrm{init}}=0.5$ and 2, respectively. In the $\beta_{i}^{\textrm{init}}=0.5$ case, the IC threshold is always at least $\sim 3$ times smaller than the mirror threshold, implying that, under realistic conditions, the IC modes continue to dominate in the $\beta_{i}^{\textrm{init}}=0.5$ case. In the $\beta_{i}^{\textrm{init}}=2$, on the other hand, the linear IC threshold is smaller than the mirror threshold by a factor $\lesssim 2$ during most of the simulations. Since we estimate that, in order for the IC instability to dominate, the mirror anisotropy threshold should be at least $\sim 2$ times larger than the one of the IC modes, this suggests that the mirror instability continues to dominate in this case.