Stochastic Ion Acceleration by the Ion-cyclotron Instability in a Growing Magnetic Field
Francisco Ley Mario Riquelme Lorenzo Sironi Daniel Verscharen Astor Sandoval

  


1Department of Astronomy, University of Wisconsin-Madison; fley@astro.wisc.edu
2Departamento de Física, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile; mario.riquelme@dfi.uchile.cl
3Department of Astronomy, Columbia University, New York, NY 10027 USA; lsironi@astro.columbia.edu
4Mullard Space Science Laboratory, University College
London, Dorking, Surrey, UK;d.verscharen@ucl.ac.uk
5Space Science Center, University of New Hampshire, Durham, NH 03824, USA
6Instituto de Astrofísica, P. Universidad Católica de Chile; asandoval@astro.puc.cl
 
 
Using 1D and 2D particle-in-cell (PIC) simulations of a plasma with a growing magnetic field B, we show that ions can be stochastically accelerated by the ion-cyclotron (IC) instability. As B grows, an ion pressure anisotropy p_⊥,i > p_||,i arises, due to the adiabatic invariance of the ion magnetic moment (p_||,i and p_⊥,i are the ion pressures parallel and perpendicular to B). When initially β_i = 0.5 (β_i ≡ 8π p_i/|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 B is amplified by a factor ∼ 2.7, this tail can be approximated as a power-law of index ∼ 3.4 plus two non-thermal bumps, and accounts for 2-3% of the ions and ∼ 18% of their kinetic energy. On the contrary, when initially β_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 β_i ≲ 1. Although we focus on cases where B is amplified by plasma shear, we check that the acceleration occurs similarly if B grows due to plasma compression. Our results are valid in a sub-relativistic regime where the ion thermal energy is ∼ 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.
 





§ 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 <cit.> and by kinetic plasma modes <cit.>. 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_⊥,i > p_||,i (where p_⊥,i and p_||,i are the ion pressures perpendicular and parallel to the local magnetic field B, respectively).

The condition p_⊥,i 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 μ_i of ions, which is defined as μ_i≡ v_⊥,i^2/B, where v_⊥,i is the ion velocity perpendicular to B and B=|B|. Thus, if B grows (decreases), the conservation of μ_i will naturally produce a pressure anisotropy with p_⊥,i > p_||,i (p_⊥,i < p_||,i). Examples of weakly collisional astrophysical plasmas where the condition p_⊥,i p_∥,i is possible are low-luminosity accretion disks around compact objects <cit.>, the intracluster medium <cit.>, and the heliosphere <cit.>.

In these systems, the growth of ion pressure anisotropy is expected to be regulated by kinetic instabilities, which break the adiabatic invariance of μ_i via pitch-angle scattering of the ions. In the p_⊥,i>p_∥,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 <cit.>. The IC instability, on the other hand, consists of propagating electromagnetic modes, with their dominant waves having k∥B <cit.>. 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_BT_i/m_ic^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 <ref> we present our simulation method and setup. In <ref> we use 2D simulations to show that ions can be accelerated by the IC instability under the condition that initially β_i ≲ 1 (β_i ≡ 8π p_i/|B|^2, where p_i is the ion isotropic pressure). In <ref> 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 <ref> 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 <ref> we show that our results are fairly independent on whether B is amplified via plasma shear or compression. In <ref> we summarize our results and present our conclusions. Additionally, in Appendix <ref> we provide details on the implementation of our 1D simulations, and in Appendix <ref> we use linear theory to analyze the applicability of our simulation results to realistic astrophysical environments.   



§ SIMULATION SETUP

 
  
llllll Parameters of the simulations  Runs   ω_c,i^init/s    m_i/m_e    β_i^init   N_ppc   L/R_L,i^init[D] 

   (or ω_c,i^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]
Simulation parameters: the initial ion cyclotron frequency ω_c,i^init (in units of s in the shearing runs and of q in the compressing runs), m_i/m_e, β_i^init, the number of particles per cell N_ppc (considering ions and electrons), and the initial box size in units of the initial ion Larmor radius L/R_L,i^init (with the number of dimensions D in squared parenthesis). In 2D, L corresponds to the height and width of the box. In all runs initially k_BT_i/m_ic^2=0.05, T_e=T_i, the electron skin depth c/ω_p,e/Δ_x=15 (where Δ_x is the grid point separation), the speed of light c=0.225 Δ_x/Δ_t (shearing runs) and 0.15Δ_x/Δ_t (compressing runs), where Δ_t is the simulation time step.  
We use the PIC code TRISTAN-MP <cit.> 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, B=B_0 x̂. This field is amplified by imposing a shear plasma velocity v = -sxŷ (represented by red arrows in Fig. <ref>a), where x is the distance along 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 d B_x/dt = dB_z/dt =0 (Fig. <ref>b shows how B changes orientation for t>0). Due to μ_j conservation, this magnetic growth drives p_⊥,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. <ref>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. <ref>b). The positions of the plasma particles are therefore given in terms of the so called `shearing coordinates', which are described both in Appendix <ref> of this paper as well as in the Appendix of <cit.>. In our 1D shearing runs, on the other hand, the simulation domain corresponds to the blue narrow stripe shown in Figs. <ref>a and <ref>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 <ref>. In <ref> 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_⊥,j/p_∥,j. For this we use the same setup as in <cit.>. In this case, B evolves as B=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 (β_i^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 (ω_c,i^init) and s (for shearing plasma runs) or q (for compressing plasma runs). The initial ion cyclotron frequency is defined as ω_c,i^init=eB_0/m_ic, 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_BT_i/m_ic^2 = 0.05. Also, these runs use m_i/m_e and ion magnetizations much smaller than expected in real astrophysical settings.[ For instance, at ∼ 10 Schwarzschild radii from the super-massive black hole Sgr A*, one expects ω_c,i/s ∼ 10^8 <cit.>, 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 ω_c,i^init/s (or ω_c,i^init/q) play a significant role.)

The numerical parameters in our runs are: the number of macro-particles per cell (N_ppc), the electron skin depth in terms of grid point spacing (c/ω_p,e/Δ_x, where ω_p,e^2=4π 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^init; R_L,i^init = v_th,i/ω_c,i^init, where v_th,i^2=k_BT_i/m_i). Table <ref> 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 <ref>.




§ 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 β_i^init=0.5), demonstrating that in this case a prominent non-thermal tail appears. Then, using simulations with β_i^init=2, we show that for β_i^init≳ 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. 


 §.§ IC vs. mirror dominated regimes
 

Figure <ref>a shows in blue-solid and red-solid lines the respective evolutions of p_⊥,i and p_∥,i for run S2m2b0.5, which uses β_i^init=0.5, m_i/m_e=2 and ω_c,i^init/s=800. We see that until t· s ≈ 1 the evolutions of both p_⊥,i and p_∥,i are in agreement with the “double adiabatic" prediction <cit.>, which is due to the conservation of μ_i and of the second adiabatic invariant. At t· s ≳ 1, the adiabatic evolution of p_⊥,i and p_∥,i is broken by the appearance of ion pressure anisotropy instabilities, which produce rapid pitch-angle scattering of the ions. 

We see from Figure <ref> that these instabilities are dominated by IC modes. Indeed, Figs. <ref>a, <ref>b and <ref>c show a snapshot at t· s=2 of the three components of δB (where δB≡B - ⟨B⟩ and ⟨⟩ denotes an average over the entire box volume). Considering that the black arrows represent the direction of ⟨B⟩, we see that δB is dominated by nearly parallel modes, which mainly appear in δ B_x and δ B_z. This is indeed consistent with the presence of transverse, circularly polarized IC modes. δ 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 δB components parallel to the plane of the simulation (they show no δ B_z component), which is in line with the expectation that δB of the mirror modes is nearly perpendicular to k×B <cit.>.



The nearly parallel IC modes can also be seen from Figs. <ref>d, <ref>e and <ref>f, which show the fluctuations in the electric field, δE.[Since our simulations are performed in the “shearing coordinate" frame <cit.>, there is not an electric field associated to the large-scale shearing motion of the plasma. Therefore, δ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_ϕ, which is related to δE by |δE| = |δB|v_ϕ/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' <cit.>, which means that their phase velocity vanishes. This implies that, as we see in Fig. <ref>, no electric field associated to the subdominant mirror modes should be present.

The dominance of the IC modes can also be seen from Figure <ref>a, which shows in solid lines the magnetic energy of δB along different axes as a function of time for run S2m2b0.5. This energy is expressed in terms of the δB components parallel to ⟨B⟩ (δ B_∥; solid-blue), perpendicular to ⟨B⟩ but parallel to the plane of the simulation (δ B_xy,⊥; solid-red), and perpendicular to both ⟨B⟩ and the plane of the simulation (δ B_z; solid-green). During most of the simulation, the energy of the magnetic fluctuations is indeed contained mainly in δ B_xy,⊥ and δ B_z, implying that the IC modes have the largest amplitude during most of the simulation time. By the end of the run (t· s=2.5), however, δ B_∥^2 becomes comparable to δ B_xy,⊥^2 and δ 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· s=2.5 (thus with a maximum B amplification factor of ∼ 2.7).

Figure <ref>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 β_i^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 β_i^init=2. We see that in this case δ B_z and δ B_xy,⊥ are subdominant in the saturated stage of the instabilities, and the energy in the magnetic fluctuations is dominated by δ B_∥^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 β_i^init∼ 1. 

In order to explore the sensitivity of this transition to m_i/m_e, in Figs. <ref>a and <ref>b we overplot δ B_∥^2, δ B_xy,⊥, and δ 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 β_i^init, k_BT_i/m_ic^2 and ω_c,i^init/s), but with m_i/m_e=10 instead of m_i/m_e=2. We see that for the two β_i^init the evolutions of δ B_xy,⊥, δ B_z, and δ B_∥ 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 β_i^init≲ 1. 

Additionally, in Appendix <ref> we use linear theory calculations to show that the condition β_i^init≲ 1 for the dominance of the IC instability should continue to hold even in realistic astrophysical plasma conditions, with m_i/m_e=1836, ω_c,i^init/s ≫ 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.
 


 §.§ 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 <cit.>

    dU_j/dt = rΔ p_j,

where r is the growth rate of the magnetic field (r ≡ (dB/dt)/B). Equation <ref> is fairly well reproduced in our simulations, as can be seen from Figure <ref>b, which shows the volume-averaged heating rate of ions dU_i/dt (solid-black) and rΔ 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 <ref>a shows the evolution of the ion spectrum for 
run S2m2b0.5 (β_i^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· s∼ 1, as seen from Figure <ref>a). By t· s=2.5 the tail can be approximated by a power law dn_i/dγ_i ∝ (γ_i-1)^-α_s with spectral index α_s ≈ 3.4 plus two bumps (γ_i is the ion Lorentz factor). The non-thermal tail at t· s=2.5 reaches Lorentz factors γ_i ∼ 10, and contains ∼ 2-3% of the ions and ∼ 18% of their energy.[After fitting the low energy part of dn_i/dγ_i to a thermal Maxwell-Boltzmann distribution, we define the non-thermal tail through the condition that dn_i/dγ_i is at least a factor 2 larger than the expectation for the thermal distribution.] The solid-black line in Fig. <ref>a represents the final (t· 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_BT_i/m_ic^2, β_i^init, and ω_c,i^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 ω_c,i^init/s will be further tested using 1D simulations in <ref>.


Figure <ref>b shows the evolution of the ion spectra for 
the 2D run S2m2b2 (m_i/m_e=2 and β_i^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· s =2.5 the tail can be approximated by a much less pronounced power law with spectral index α_s ≈ 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 <ref> we show that this is indeed the case making use of 1D simulations in which the mirror modes are artificially suppressed. In <ref> 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 ω_c,i^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.



§ 1D SHEARING SIMULATIONS
 

Since the IC modes propagate mainly parallel to the background magnetic field ⟨B⟩, in this section we study the ion acceleration due to the IC instability by only capturing modes with wave vector k parallel to ⟨B⟩. We do so by using the 1D version of our shear coordinates simulations <cit.>. 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 ⟨B⟩. A detailed description of our 1D setup is in Appendix <ref>. 

In the next section, we show the suitability of the 1D setup to study problems dominated by modes with k∥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 ω_c,i^init/s, which will make use of the low computational cost of the 1D runs.



 §.§ 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. <ref>a shows the evolution of the ion spectrum for run S1m2b0.5 (1D) from t· s=0 to t· 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 <ref>a. Indeed, both spectra can be described as a power law of index α_s≈ 3.4 plus two bumps, with the intermediate energy bump appearing at t· s ≳ 2. A more detailed comparison can be made by overplotting the final (t· s=2.5) spectrum for run S2m2b0.5 in Figure Fig. <ref>a (black line). We see that the two final spectra are very similar, with the main difference being a ∼ 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. <ref>a shows p_i,∥ (red) and p_i,⊥ (blue) for runs S1m2b0.5 (1D; dotted) and S2m2b0.5 (2D; solid). Both for p_i,∥ and p_i,⊥, the 1D and 2D simulations give essentially the same results. In dashed lines we show the corresponding double-adiabatic behavior <cit.>, which is followed quite well by the two  simulations until t· s ≈ 1. With respect to energy conservation, Figure <ref>b shows the volume average dU_i/dt (black) and rΔ 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 δB, we use δ 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 <ref>a shows a snapshot of δ B_z at t· s=2 for the 1D run S1m2b0.5. We see that both in terms of the dominant wavelength (≈ 30 R_L,i^init) and of its amplitude, δ B_z behaves fairly similarly to the 2D case, shown in Fig. <ref>c. Figures <ref>b, <ref>c and <ref>d show power spectra of δ 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· s < 1.2, 1.8< t· s < 2, and 2.2< t· 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 ∼ 1.5) shift in the peak of the 1D spectra towards longer wavelengths, and a factor ∼ 2 increase in the peak amplitude.

These differences in the δ B_z spectra can be explained to a large extent by the small differences in the ion energy spectra. Indeed, in <ref> 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 ∝ 1/γ_i. Thus, since the 1D runs produce maximum ion energies ∼ 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 δ B_z. For relativistic ions interacting resonantly with parallelly propagating waves, the effective scattering frequency ν_eff,i should scale as <cit.>:

    ν_eff,i∝ω_c,i/B_0^2γ_id(δ B_z^2)/dln(k).
        
 
Thus, since the ions in the 1D and 2D runs are scattered at roughly the same rate (given their similar evolution of p_∥,i and p_⊥,i, as shown in Fig. <ref>a), Eq. <ref> implies that the peak value of dB_z^2k/dk in the 1D case should roughly be ∼ 2 times larger than in the 2D case, which is seen in panels b, c and d of Fig. <ref>.

Thus, besides a factor ∼ 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 ω_c,i^init/s.



 §.§ 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 (β_i^init=2) with an analogous 1D run where these modes are artificially suppressed. This is done in Fig. <ref>b, which shows in black the final spectrum of the 2D run S2m2b2, where the mirror modes dominate. In addition, Fig. <ref>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. <ref>. Thus, when the mirror modes dominate, the scattering of ions tends to be elastic, and no acceleration effect should be present.



 §.§ m_i/m_e and ω_c,i^init dependence




We use 1D simulations to explore the dependence of the ion acceleration by the IC instability on m_i/m_e and ω_c,i^init/s, focusing on the case with β_i^init=0.5. In terms of the evolutions of p_⊥,i and p_∥,i, Figure <ref>a shows the cases of runs S1m8b0.5 (ω_c,i^init/s=800, m_i/m_e=8; dotted line), S1m128b0.5 (ω_c,i^init/s=800, m_i/m_e=128; solid line) and S1m8b0.5d (ω_c,i^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 <ref>b shows dU_i/dt (black) and rΔ p_i (green) for the same runs. We see that the ion energy gain is fairly independent of m_i/m_e and ω_c,i^init/s, and in all cases it reasonably well agrees with the heating prediction through anisotropic viscosity.

Figure <ref>a shows the final ion spectra (t· 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 <ref>). In all simulations the ions share the same parameters: ω_c,i^init/s=800, β_i^init=0.5, and k_BT_i=k_BT_e=0.05m_ic^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 α_s ≈ 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 <cit.>. 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 <ref>b shows the final spectra for simulations with m_i/m_e=8, β_i^init=0.5, and k_BT_i=k_BT_e=0.05m_ic^2, but with ω_c,i^init/s=400, 800, 1600 and 3200 (runs S1m8b0.5b, S1m8b0.5, S1m8b0.5c and S1m8b0.5d, respectively). The spectra get slightly harder as ω_c,i^init/s increases, with the difference between them being progressively less significant as ω_c,i^init/s grows. However, in all cases the tail can be well described as a power-law of spectral index α_s ≈ 3.4 plus two bumps. 

The independence of the ion acceleration on ω_c,i^init/s can be inferred from the way the effective ion scattering rate ν_eff,i provided by the IC waves is related to s. This scattering rate can be estimated from the evolution of p_∥,i in a homogeneous plasma with no heat flux, assuming that B evolves on time and length scales much larger than ω_c,i^-1 and R_L,i (which is the case in our runs). This evolution is given by Eq. 1 of <cit.>:

    ∂ p_∥,i/∂ t + ∇· (p_∥,iv) + 2p_∥,ib̂b̂:∇v=2/3ν_eff,iΔ p,

where v is the plasma bulk velocity and b̂≡B/B. In the case of the shearing plasma (v=-sxŷ), ∇· (p_∥,iv)=0 and b̂b̂:∇v=b̂_xb̂_ys. Fig. <ref>a shows that, after the saturation of the IC modes, p_∥,i changes at a rate much smaller than s so we can approximate ∂ p_∥,i/∂ t ≈ 0. Additionally, b̂_xb̂_y ranges between 0.5 at t· s=1 and 0.34 at t· s=2.5. Thus, we simply assume b̂_xb̂_y≈ 1/2 and Eq. <ref> becomes:

    ν_eff,i≈3/2s (p_∥,i/Δ p_i).

By comparing the evolutions of p_⊥,i and p_∥,i for runs S1m8b0.8 (ω_c,i^init/s=800; dot-dashed line) and S1m8b0.8d (ω_c,i^init/s=3200; dashed line), Figure <ref>a shows that the factor p_∥,i/Δ p_i is fairly independent of ω_c,i^init/s. Thus, Eq. <ref> implies that ν_eff,i∝ s.

This proportionality between ν_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. <ref>b), should make the accelerating effect of the IC modes independent of ω_c,i^init/s.

Our 1D simulations, therefore, show that the ion acceleration by IC modes is fairly independent of m_i/m_e and ω_c,i^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.



§ THE ACCELERATION MECHANISM
 

The two possible sources of energy for the ions in our simulations are: the energy gain due to anisotropic viscosity (Eq. <ref>), 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 <ref> we identify the contributions of each of these energy sources to producing the non-thermal ion spectra.



 §.§ IC work vs. anisotropic viscosity
 

Figure <ref> 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· 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. <ref>a. These ions are chosen so that at t· s=2.5 their Lorentz factors satisfy 0.09 < γ_i-1 < 0.1 (marked by the vertical red line in Fig. <ref>a).   

  * The ions in the “intermediate energy" bump of the tail (shown in Fig. <ref>b). These ions have Lorentz factors in the range 1 < γ_i-1 < 1.03 at t· s=2.5 (marked by the vertical dark green line in Fig. <ref>a).

  * The ions in the “high energy" bump of the tail (shown in Fig. <ref>c). This population corresponds to the highest energy ions, defined by γ_i-1 > 4.8 at t· s=2.5 (marked by the grey region in Fig. <ref>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. <ref>.

  * The energy gain by anisotropic viscosity, AV, shown by the green line in Fig. <ref>.[For each ion population, this energy gain is calculated as the integral in time of the rate of energy gain due to viscosity: ∫ dt  r Δ p_i (see Eq. <ref>), where Δ p_i is calculated using only the ions of each population.] 


The blue line in Figure <ref>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 <ref>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· s ≈ 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 <ref>a. Finally, Figure <ref>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 <ref>a, <ref>b and <ref>c we also plot the sum of AV and W_IC (red line) and the average change in energy of the three ion populations, ⟨Δ E_i ⟩ (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.   



 §.§ Resonance with IC waves
 


In <ref> 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:

    ω/k - v_|| = ω_c,i/γ_ik,

where ω 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 ω_c,i (≡ eB/m_ic) 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. <ref>a and <ref>b we measure the ratio ω/k for the IC waves at two different time intervals for the 1D run S1m2b0.5 (β_i^init=0.5, m_i/m_e=2, ω_c,i^init=800). We do this by plotting |δB̃_z(ω,k) + i δB̃_⊥, xy(ω,k)|^2, where δ B_z(ω,k) and δ B_⊥, xy correspond to two mutually perpendicular components of δB, that are also perpendicular to ⟨B⟩ (see caption of Fig. <ref>), the tilde () denotes the time and space Fourier transform of a quantity and i=√(-1). Thus, the combination δB̃_z(ω,k) + i δB̃_⊥, xy(ω,k) allows to separate the contributions to δ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 <cit.>. In the case of run S1m2b0.5, the IC waves contribute to |δB̃_z(ω,k) + i δB̃_⊥, xy(ω,k)| only for ω > 0, while the whistler wave contribution appears for ω < 0. This way, calculating |δB̃_z(ω,k) + i δB̃_⊥, xy(ω,k)|^2 allows to separate the IC and whistler contributions to δB, and to estimate ω(k) for these two modes.

Figure <ref>a corresponds to the time interval 1.1 < t· s < 1.2 of run S1m2b0.5. The IC modes have a phase velocity of ω/k ≈ 0.5c. Additionally, we measure the rms ion velocity parallel to B at 1.1 < t· s < 1.2, which is v_∥^rms≈ 0.16c,[This can be estimated by the factor ∼ 2 decrease in p_i,∥ seen in Fig. <ref>a, and considering that initially k_BT_i/m_ic^2=0.05.] implying that to a good approximation we can neglect the v_∥ term on the left hand side of Eq. <ref>. The resonance condition at 1.1 < t· s < 1.2 can thus be written as:


    γ_i(kc/ω_c,i) ≈ 2.
 
At 1.1 < t· s < 1.2, γ_i is in the range 1 < γ_i≲ 2 (see Fig. <ref>a), implying that most of the power of the IC modes should be in the range 1≲ kc/ω_c,i≲ 2, which coincides well with the range of k in which most of the IC power is observed in Fig. <ref>a. 

Analogously, Fig. <ref>b shows |δ B_z(ω,k) + i δ B_⊥, xy(ω,k)|^2 for the same simulation but in the time range 2.2 < t· s < 2.3. In this case, Fig. <ref>b shows that ω/k ≈ 0.7c. Thus, making a similar analysis as in the case 1.1 < t· s < 1.2, we obtain:[Here we also assume ω/k (≈ 0.7c) ≫ v_∥, which allows neglecting v_∥ in Eq. <ref>. This is a reasonable approximation considering that at 2.2 < t· s < 2.3, v_∥ (measured directly from the simulation) is always smaller than ∼ 0.4c, even considering the highest energy particles in the non-thermal tail.]

    γ_i(k c/ω_c,i) ≈ 1.5.

Since in this time interval 1 < γ_i≲ 8, we obtain that most of the power of the IC modes should be in the range 0.2 ≲ kc/ω_c,i≲ 1.5. This interval coincides reasonably well with the range of k where most power is concentrated in Fig. <ref>b. Notice that this power appears to be enhanced in two intervals of k. The high-k interval corresponds to 0.7≲ kc/ω_c,i≲ 1.5, which, according to Eq. <ref>, implies resonance with ions with 1 < γ_i≲ 2. Remarkably, this is a range of γ_i with abundant IC scattering at 2.2 ≲ t· s ≲ 2.3, as shown by the rapid formation of the “intermediate energy" bump, which mainly occurs at 2 ≲ t· s ≲ 2.5. The low-k enhancement occurs for 0.2≲ kc/ω_c,i≲ 0.5, which, according to Eq. <ref>, corresponds to 3 ≲γ_i ≲ 8. This γ_i interval coincides well with the “high energy" bump shown in Fig. <ref>a at t· s ∼ 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 <ref>c and <ref>d show the same quantities as Figs <ref>a and <ref>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 |δB̃_z(ω,k) + i δB̃_⊥, xy(ω,k)|^2 essentially preserves the same properties for the IC modes (ω > 0). The fact that ω(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 ω_c,e/ω_p,e, ω_c,e/s and k_BT_e/m_ec^2 do change when varying m_i/m_e.

Finally, the four panels in Fig. <ref> also include IC and whistler theoretical dispersion relation calculations, ω_theo(k) (dashed lines), obtained with the linear Vlasov solver NHDS <cit.>. 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 ω_theo(k) reproduces the behavior of |δB̃_z(ω,k) + i δB̃_⊥, xy(ω,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 γ_e ≫ 1), which is a regime strictly not captured by NHDS.



§ 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 ⟨B⟩, 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 ⟨B⟩ does not play an important role in the acceleration. We do this by running 1D simulations in which ⟨B⟩ grows due to plasma compression. The simulation setup is the same as in <cit.>. The plasma is compressed along the y and z axes, with ⟨B⟩ pointing along x. In this setup, ⟨B⟩ is amplified due to magnetic flux conservation <cit.>. In this situation, the background field grows as |⟨B⟩|=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 k∥⟨B⟩, 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: β_i^init=0.5, k_BT_i/m_ic^2=0.05, and T_i=T_e, which allows a direct comparison between the two setups. Fig. <ref>a shows the evolution of T_⊥,i (solid-blue) and T_∥,i (solid-red) for run C1m8b0.5a, which uses m_i/m_e=8 and ω_c,i^init/q=1600 (k_BT_⊥,i=p_⊥,i/n_i and k_BT_∥,i=p_∥,i/n_i, where n_i is the ion density).[We show T_⊥,i and T_∥,i instead of p_⊥,i and p_∥,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 <cit.>.[The small initial discrepancies (t· s ≲ 0.6) between the double-adiabatic predictions and T_∥,i and T_⊥,i are because these predictions assume that the ions are non-relativistic, which is mildly broken for our k_BT_i/m_ic^2 ∼ 0.1 ions.] The adiabatic evolution breaks at t· q ≈ 0.6, after the exponential growth of the IC modes begins. This exponential growth and subsequent saturation can be seen from Fig. <ref>b, which shows the δ B_y^2 and δ B_z^2 components of run C1m8b0.5a, which evolve quite similarly to the δ B_⊥,xy and δ B_z^2 fluctuations of analogous shearing runs (see, e.g., the cases of runs S2m2b0.5 and S2m10b0.5 in Fig. <ref>a).

The evolution of the ion spectrum of run C1m8b0.5a is shown in Fig. <ref>c. We see a remarkable similarity with the spectral evolution of the shearing runs with β_i^init=0.5, shown in Figs. <ref>a and <ref>a. Indeed, also in the compressing case, there is a growth of a non-thermal tail that, by t· q=2.5, can be described as a power law of index α≈ 3.5 plus two bumps. 

We also test the dependence of this acceleration on both m_i/m_e and ion magnetization ω_c,i^init/q. Fig. <ref>d shows a comparison between cases with m_i/m_e=8 and 16 (black and red, respectively), and with ω_c,i/q=1600 and 3200 (solid and dashed, respectively), which share the same initial parameters: β_i^init=0.5, k_BT_i/m_ic^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 ω_c,i/q increases, in agreement with the shearing results of <ref>. 

These results essentially reproduce our shearing runs, both in terms of the significance of the ion acceleration in the regime: β_i^init=0.5, k_BT_i/m_ic^2=0.05 and T_i=T_e, and in terms of the almost independence on m_i/m_e and ω_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 ω_c,i^init/q implies, similarly to the shearing case, that the effective ion scattering rate ν_eff,i has to be ∝ q. This can indeed be inferred by applying Eq. <ref> to the compressing runs. In this case, b̂b̂:∇v=0, thus (using the continuity relation ∂ n_i/∂ t = -n_i∇·v) n_i∂ T_∥,i/∂ t = (2/3) ν_eff,iΔ p_i. However, Fig. <ref>a shows that, after the saturation of the IC waves, ∂ T_∥,i/∂ t ∼ (q/2) T_∥,i, which suggests:

    ν_eff,i≈3/4q (p_∥,i/Δ p_i).

The ratio p_∥,i/Δ p_i evolves fairly independent of ω_c,i^init/q, as can be seen from Fig. <ref>a, which shows T_⊥,i and T_∥,i for runs C1m8b0.5a (ω_c,i^init/q=1600; solid line) and C1m8b0.5b (ω_c,i^init/q=3200; dot-dashed line). Therefore, Eq. <ref> implies that ν_eff,i∝ q.

The similarity between Eqs. <ref> and <ref>, 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_∥,i/Δ p_i∼ 1/2 during most of the saturated IC regime (see Figs. <ref>a and <ref>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.[Notice, 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 β_i^init=0.5 the final value of the parallel ion beta is β_∥,i^final∼ 0.025 (considering the evolution of p_∥,i seen in Fig. <ref>a and the expected evolution of B), while in the analogous compressing runs β_∥,i^final∼ 0.1 (considering the evolution of T_∥,i seen in Fig. <ref>a and the expected evolutions of n_i and B). This, plus the different factors on the right hand sides of Eqs. <ref> and <ref>, imply that the final ion spectra in these two setups should be similar but not necessarily the same.]



§ 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 ⟨B⟩.

In the regime in which initially k_BT_i=k_BT_e=0.05m_ic^2, we show that the IC ion acceleration is significant if β_i^init≲ 1. This is demonstrated by comparing 2D simulations with β_i^init=0.5 and 2. When β_i^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 ∼ 2.7, the tail can be approximately described as a power-law of index ∼ 3.4 plus two non-thermal bumps. The maximum ion Lorentz factor at that time is γ_i ∼ 10, but it continues to grow at the end of the simulation. Also, the tail accounts for 2-3% of the ions and ∼ 18% of their kinetic energy. On the other hand, when β_i^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. <ref>). 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 <cit.>. 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 <ref>). 

Given that our simulations can not use realistic values of m_i/m_e and ω_c,i^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 β_i^init=0.5 does not depend on m_i/m_e and ω_c,i^init/s. This was done comparing 2D simulations with m_i/m_e=2 and 10 (see <ref>), 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 <ref>). Both analyses show that, in realistic astrophysical plasmas, the IC instability dominates in the regime β_i^init≲ 1, at least for the case of T_e=T_i explored here. 

Then, using that the dominant IC wave vectors k satisfy k∥⟨B⟩, 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 ω_c,i^init/s. While m_i/m_e almost does not affect the acceleration, increasing ω_c,i^init/s only produces a slight hardening of the non-thermal tail (see <ref>). This almost complete independence of the acceleration on ω_c,i^init/s is consistent with the effective ion scattering rate ν_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 <cit.>. We find essentially no difference in the ion acceleration efficiency between the shearing and compressing cases (see <ref>).

Our work is valid in a sub-relativistic regime in which initially k_BT_i=k_BT_e=0.05m_ic^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 β_i^init≲ 1 required for the acceleration is most likely satisfied in the coronal region of the disks <cit.>. 

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_BT_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 ≪ 1 may increase significantly the values of β_i^init for which the IC instability dominates <cit.>. Since the condition T_e/T_i ≪ 1 is most likely satisfied in low luminosity disks <cit.>, 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 δ B/B∼ 1 amplifications and decreases of the field <cit.>. 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.

99
tocsectionBibliography

[Anderson et al. (1991)]AndersonEtAl1991 Anderson, B. J., Fuselier, S. A., & Murr, D. 1991, GeoRL, 18, 1955
[Bale et al. (2009)]BaleEtAl2009 Bale, S. D.,  Kasper, J. C., Howes, G. G., Quataert, E., Salem, C., & Sundkvist, D. 2009, PRL 103, 211101
[Buneman(1993)]Buneman93 Buneman, O. 1993, “Computer Space Plasma Physics”, Terra Scientific, Tokyo, 67
[Chael et al. (2018)]ChaelEtAl2018 Chael, A, Rowan, M, Narayan, R., Johnson, M. & Sironi, L. 2018, MNRAS, 478, 5209
[Chandran (2003)]Chandran2003 Chandran, B. 2003, ,  599, 1426.
[Chew et al. (1956)]ChewEtAl1956 Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, RSPSA, 236, 112
[Cho & Lazarian (2006)]ChoEtAl2006 Cho, J. & Lazarian, A. 2006, , 638, 811
[Dermer et al. (1996)]DermerEtAl1996 Dermer, C. D., Miller, J. A. & Li, H. 1996, , 456, 106
[Gary (1992)]Gary1992 Gary, S. P. 1992, JGR, 97, 8519
[Gary & Wang (1996)]GaryEtAl1996 Gary, S. P. & Wang, J. 1996, J. Geophys. Res., 101, 10749
[Hasegawa (1969)]Hasegawa1969 Hasegawa, A. 1969, PhFl, 12, 2642
[Hellinger et al. (2006)]HellingerEtAl2006 Hellinger, P., Travnicek, P.,  Kasper, J. C., &  Lazarus, A. J. 2006, GRL, 33, L09101
[Isenberg et al. (2013)]IsenbergEtAl2013 Isenberg, P. A., Maruca, B. A., & Kasper, J. C. 2013, ApJ, 773, 164
[Kulsrud & Pearce (1969)]KulsrudEtAl1969 Kulsrud, R. M. & Pearce, W. P. 1969, , 156, 445
[Kulsrud(1983)]Kulsrud1983 Kulsrud, R. M. 1983, in Handbook of Plasma Physics, ed. M. N. Rosenbluth & R. Z. Sagdeev (Amsterdam: North Holland), 115
[Lynn et al.(2014)]LynnEtAl2014 Lynn, J., Quataert, E., & Chandran, B., & Parrish, I. 2014, 791, 71
[Lyutikov (2007)]Lyutikov2007 Lyutikov, M. 2007, ApJL, 668, L1
[Maruca et al.(2011)]MarucaEtAl2011 Maruca, B. A., Kasper, J. C., & Bale, S. D. 2011, PhRvL, 107, 201101
[Narayan & Yi (1995)]NarayanEtAl1995 Narayan, R. & Yi, I. 1995, , 452, 710
[Petrosian & Liu (2004)]PetrosianEtAl2004 Petrosian, V. & Liu, S. 2004, , 610, 550
[Pokhotelov et al.(2004)]PokhotelovEtAl2004 Pokhotelov, O. A., Sagdeev,R. Z., Balikhin, M. A.,& Treuman,R.A. 2004, JGR, 109, A09213
[Ponti et al.(2017)]PontiEtAl2017 Ponti, G.,  George, E.,  Scaringi, S., Zhang, S., Jin, C., et al. 2017, , 468, 2447
[Riquelme et al.(2012)]RiquelmeEtAl2012 Riquelme, M. A., Quataert, E., Sharma, P., & Spitkovsky, A. 2012, , 755, 50 
[Riquelme et al.(2017)]RiquelmeEtAl2017 Riquelme, M. A., Osorio, A., & Quataert, E. 2017, , 850, 113
[Schekochihin et al.(2005)]SchekochihinEtAl2005 Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Hammett, G. W., & Sharma, P. 2005, ApJ, 629, 139
[Sharma et al.(2006)]SharmaEtAl2006 Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. 2006, ApJ, 637, 952
[Sharma et al.(2007)]SharmaEtAl2007 Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. 2007, , 667, 714
[Sironi & Narayan (2015)]SironiEtAl2015 Sironi, L., & Narayan, R. 2015, ApJ, 800, 88
[Sironi (2015)]Sironi2015 Sironi, L. 2015, ApJ, 800, 89
[Snyder et al.(1997)]SnyderEtAl1997 Snyder, P. B., Hammett, G. W., & Dorland, W. 1997, Phys. Plasmas, 4, 3974 
[Southwood & Kivelson (1993)]SouthwoodEtAl1993 Southwood, D. J., & Kivelson, M. G. 1993, JGR, 98, 9181
[Spitkovsky(2005)]Spitkovsky05 Spitkovsky, A. 2005, AIP Conf. Proc, 801, 345, astro-ph/0603211
[Verscharen at al.(2016)]VerscharenEtAl2016 Verscharen, D., Chandran, B., Klein, K. G., & Quataert, E. 2016, 831, 128
[Verscharen & Chandran (2018)]VerscharenEtAl2018 Verscharen, D. & Chandran, B. 2018, Res. Notes AAS 2, 13
[Verscharen et al. (2018)]VerscharenEtAl2018b Verscharen, D., Klein, K. G., Chandran, B., Stevens, M. L., Salem, C. S., & Bale, S. D. 2018, J. Plasma Phys., 84, 905840403
[Verscharen et al. (2019)]VerscharenEtAl2019 Verscharen, D., Klein, K. G., & Maruca, B. A. 2019, arXiv:1902.03448
[Yuan et al.(2003)]YuanEtAl2003 Yuan, F., Quataert, E., & Narayan, R. 2003, , 598, 301




§ 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 ⟨B⟩. This is equivalent to assuming that the wave vectors k captured in the simulations satisfy k∥⟨B⟩. Since in a shearing plasma the direction of ⟨B⟩ 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 <cit.>. 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ŷ). Figures <ref>a and <ref>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',y',z',t') are defined in terms of the regular, inertial coordinates (x, y, z, t) as 

    x'=x,  y'=Γ(y-vt),  z'=z, and t'=Γ(t-vy/c^2),

where Γ=(1-v^2/c^2)^-1/2, v=-c arctanh(sx/c),[Notice that, in the limit |v| ≪ c, the v=-c arctanh(sx/c) expression is equivalent to v=-sx, as we assume in Figures <ref>a and <ref>b. Since v=-c arctanh(sx/c) ensures that |v| < c, this expression was adopted in <cit.> in order to deduce the equations that describe the dynamics of the plasma in the shearing coordinate system. However, both in <cit.> and in this paper we are interested in the plasma dynamics in the regime |v| ≪ c.] c is the speed of light and s is the shear rate of the plasma. 

Our 1D domain is defined by y'=0, and it is shown by the blue lines in Figures <ref>a and <ref>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 ⟨B⟩ will always be parallel to the 1D domain at t>0, as also depicted in Figures <ref>a and <ref>b. Therefore, if the problem of interest is dominated by waves that propagate parallel to ⟨B⟩, our 1D simulations will be able to capture the essence of the phenomenon.

In order for our 1D simulations to consistently satisfy k∥⟨B⟩, we need to replace x' with a new coordinate x'_1, which we define as 

    x'_1≡ x'-y'st'/1+s^2t'^2.

Indeed, if we assume |v| ≪ c and sy' ≪ c, the partial derivatives with respect to the spatial coordinates in the inertial frame are

    ∂/∂ x = 1/1+s^2t^2∂/∂ x'_1 + st∂/∂ y', ∂/∂ y = -st/1+s^2t^2∂/∂ x'_1 + ∂/∂ y', and∂/∂ z = ∂/∂ z'.

Thus, if we impose the 1D condition, namely that the fields depend only on x'_1, with ∂/∂ y' = ∂/∂ z' = 0, Equations <ref> 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 ⟨B⟩ 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 ∂/∂ y' = ∂/∂ z' = 0, our 1D simulations are able to capture all the modes with wave vectors k parallel to ⟨B⟩, despite the fact that ⟨B⟩ changes orientation over time. 

The assumptions |v| ≪ c and sy' ≪ c are equivalent to assuming that the plasma region of interest has a typical size, L, that satisfies L ≪ 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 ≪ c/s applied to species j becomes: [γ_js/ω_c,j][v_j,⊥/c]≪ 1, where γ_j, ω_c,j, and v_j,⊥ are the Lorentz factor, (non-relativistic) cyclotron frequency, and velocity perpendicular to B of particles j, respectively. In typical astrophysical environments, ω_c,j is many orders of magnitude larger than s, therefore the assumption L ≪ 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 ≪ c/s is assumed in all of these expressions.     



 §.§ Evolution of E' and B'


Considering the evolutions of E' and B' in terms of J' in the shearing coordinates presented in Equations A14 and A26 of <cit.>, after replacing x' with x'_1 (Eq. <ref>) and assuming ∂/∂ y' = ∂/∂ z' = 0, the fields dynamics is given by:

    ∂ B_x'/∂ t' = cst'/1+s^2t'^2∂ E_z'/∂ x'_1, ∂ B_y'/∂ t' = c1/1+s^2t'^2∂ E_z'/∂ x'_1 - sB_x', ∂ B_z'/∂ t' = -c/1+s^2t'^2( ∂ E_y'/∂ x'_1 + st'∂ E_x'/∂ x'_1),

and

    ∂ E_x'/∂ t' = -cst'/1+s^2t'^2∂ B_z'/∂ x'_1-4π J'_x, ∂ E_y'/∂ t' = -c1/1+s^2t'^2∂ B_z'/∂ x'_1 - sE_x'-4π J'_y, ∂ E_z'/∂ t' = c/1+s^2t'^2(∂ B_y'/∂ x'_1 + st'∂ B_x'/∂ x'_1) -4π J'_z.

In Equations <ref> and <ref>, E', B', and J', 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 <cit.>). This means that J' can be calculated directly from the motion of the particles in the shear coordinates, assuming that their momenta p' can be obtained from their momenta in the inertial frame p, using the conventional relativistic momentum transformation (see Eq. A30 of <cit.>). 



 §.§ Evolution of particle momenta

Since our 1D setup uses the same definitions of time (t') and of particle momenta (p') as in the shearing coordinate system of <cit.>, the evolution of p' will not change and will be given by

    dp'/dt' = q_c(E' + u'/c×B') + sp_x'ŷ,

where q_c is the particles' electric charge and u'=p'/mγ', with m and γ' being the particle mass and Lorentz factor, respectively. Eq. <ref> corresponds to Eqs. A31 of <cit.> in the limit L ≪ c/s. 



 §.§ Evolution of particle positions
  
By taking the derivative of Equation <ref> with respect to t', and applying the limit L ≪ c/s to the expressions for dx'/dt' and dy'/dt' given by Eq. A35 of <cit.>, we find  

    dx'_1/dt' = u_x'1/1+s^2t'^2 - u_y'st'/1+s^2t'^2,
 
where u_x' and u_y' are the x and y components of u'. Eq. <ref> implies that dx'_1/dt' is the scalar product between u' and ⟨B⟩/|⟨B⟩| (=(x̂ - st'ŷ)/(1+s^2t'^2)^1/2), corrected by the “expansion" of the 1D domain (represented by the y'=0 region depicted in blue in Figures <ref>a and <ref>b), which provides the extra common factor 1/(1+s^2t'^2)^1/2 on the right hand side of Eq. <ref>. This implies that x'_1 correctly represents particle displacements along ⟨B⟩.



 §.§ 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: ∇·E=4πρ_c, where E and ρ_c are the electric field and charge density in that frame. Indeed, applying the transformations given by Eqs. <ref> to Gauss' law in the limit L ≪ c/s, and assuming first the more general problem in 3D described in <cit.>:

    ∂ρ_c'/∂ t' + ∂ J_x'/∂ x' + ∂ (J_y'+st'J_x') /∂ y' + ∂ J_z'/∂ z'=0,

where ρ_c' is the charge density of the plasma calculated using the shearing coordinates, which, as shown in Eq. A29 of <cit.>, is equal to ρ_c in the limit L ≪ c/s. Thus, according to Eq. <ref>, in order to satisfy charge conservation in our multidimensional shearing coordinate, particle motions in y' lead to an effective current along the y' direction of J'_y+st'J'_x. Thus, charge conservation is ensured in our multidimensional shearing coordinate simulations by evolving x', y' and z' according to


    dx'/dt'=u_x', dy'/dt'=u_y'+st'u_x', anddz'/dt'=u_z',

which we do in the 2D simulations presented in this paper and in our previous works <cit.>. Eq. <ref> is equivalent to Eq. A35 of <cit.> in the limit L ≪ c/s.[Notice that the notation used in Eq. A35 of <cit.> is different from the one used in this work, and that we interprete their u_x', u_y' and u_z' as our dx'/dt', dy'/dt' and dz'/dt'.]

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

    ∂ρ_c'/∂ t' + ∂/∂ x'_1(J_x'/1+s^2t'^2 - st'J_y'/1+s^2t'^2)=0.

Thus, according to Eq. <ref>, in order for our 1D runs to conserve charge, x'_1 must evolve so that the effective current along ⟨B⟩ is (J_x'-st'J_y')/(1+s^2t'^2), which is indeed satisfied by the evolution of x_1' given by Eq. <ref>. 



§ B. IC VS. MIRROR UNDER REALISTIC CONDITIONS: LINEAR THEORY ANALYSIS



In order to inquire whether the dominance of the IC instability for β_i^init≲ 1 continues to be valid in realistic, astrophysically relevant regimes (m_i/m_e=1836 and ω_c,i/s ≫ 800, which we can not study with our 2D simulations), we make use of linear theory. Thus, we calculate the threshold ion pressure anisotropy, Δ p_i/p_∥,i≡ (p_⊥,i-p_∥,i)/p_∥,i, needed for the growth of IC and mirror modes at a rate s, using the linear Vlasov solver NHDS <cit.>.[Modern linear solvers to the linear Vlasov-Maxwell system of equations can account for relativistic effects and non-Maxwellian background distributions <cit.>. 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 Δ p_i/p_∥,i at a nearly stationary level (as we see in Figs. <ref>a and <ref>b), the modes that provide the pitch-angle scattering must grow at roughly the rate at which the anisotropy is driven, which is ∼ 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 <cit.>. 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 <cit.>. Another possibility is the inhomogeneity in the magnetic field produced by the (subdominant but still present) mirror modes, which may also affect this resonance <cit.>. 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 <ref>a and <ref>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 β_i^init=0.5 and 2, respectively (both with m_i/m_e=2 and ω_c,i^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 ∼ 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 ∼ 2 times smaller than the mirror threshold (Fig. <ref>a), while in the mirror dominated case (Fig. <ref>b) the IC threshold is at most ∼ 1.5 times smaller. These results suggest that, in order for the IC instability to dominate, the IC threshold should be at least ∼ 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ω_c,i. Figures <ref>a and <ref>b show in dotted lines the corresponding linear anisotropy thresholds for the mirror (dotted red) and IC (dotted blue) instabilities for β_i^init=0.5 and 2, respectively. In the β_i^init=0.5 case, the IC threshold is always at least ∼ 3 times smaller than the mirror threshold, implying that, under realistic conditions, the IC modes continue to dominate in the β_i^init=0.5 case. In the β_i^init=2, on the other hand, the linear IC threshold is smaller than the mirror threshold by a factor ≲ 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 ∼ 2 times larger than the one of the IC modes, this suggests that the mirror instability continues to dominate in this case.