PIC Simulations of Continuously Driven Mirror and Ion Cyclotron Instabilities in High Beta Astrophysical and Heliospheric Plasmas

1.1 MB·

by  Mario Riquelme

Draftversion January13,2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
PIC SIMULATIONS OF CONTINUOUSLY DRIVEN MIRROR AND ION CYCLOTRON INSTABILITIES IN
HIGH BETA ASTROPHYSICAL AND HELIOSPHERIC PLASMAS
Mario A. Riquelme1, Eliot Quataert2 & Daniel Verscharen3
Draft version January 13, 2015
ABSTRACT
We use particle-in-cell (PIC) simulations to study the nonlinear evolution of ion velocity space
5
instabilities in an idealized problem in which a background velocity shear continuously amplifies
1
the magnetic field. We simulate the astrophysically relevant regime where the shear timescale is
0
long compared to the ion cyclotron period, and the plasma beta is β 1 100. The background
2 ∼ −
field amplification in our calculation is meant to mimic processes such as turbulent fluctuations

n or MHD-scale instabilities. The field amplification continuously drives a pressure anisotropy with
a p > p and the plasma becomes unstable to the mirror and ion cyclotron instabilities. In all cases,
⊥ k
J the nonlinear state is dominated by the mirror instability, not the ion cyclotron instability, and the

0 plasmapressureanisotropysaturatesnearthethresholdforthelinearmirrorinstability. Themagnetic
1 field fluctuations initially undergo exponential growth but saturate in a secular phase in which the
fluctuations grow on the same timescale as the background magnetic field (with δB 0.3 B in the

] secular phase). At early times, the ion magnetic moment is well-conserved but once∼the flhucituation
E
amplitudes exceedδB 0.1 B , the magnetic momentis no longerconservedbut insteadchangeson
H ∼ h i
a timescale comparable to that of the meanmagnetic field. We discuss the implications of ourresults
. for low-collisionality astrophysical plasmas, including the near-Earth solar wind and low-luminosity
h
accretion disks around black holes.
p
-
o Subject headings: accretion, accretion disks – instabilities – plasmas – solar wind
r
t
s
a 1. INTRODUCTION Velocity space instabilities limit the growth of the pres-
[ sure anisotropy and give rise to small scale fluctuations
Ion pressure anisotropies are ubiquitous in heliospheric
2 and astrophysical plasmas. In the absence of Coulomb in the magnetic field. These fluctuations can, in turn,
v collisions, the magnetic moment of ions, µi (≡ v⊥2,i/B, abffyemctotdhifeyilnagrgtehsecamleeatnrafnrespeoprattphroofpepratriteiscloefs.the plasma
4 where v is the ion velocity perpendicular to the
⊥,i
1 local magnetic field B~ and B = B~ ) is an adiabatic Magnetic field amplification in a collisionless system
0 | |
invariant. Thus, if ion collisions are infrequent, the ion generically drives p > p . There are two ion-scale
0 ⊥,i ||,i
velocity distributions parallel and perpendicular to the instabilities that can be excited in this regime: the mir-
.
2 magnetic field decouple, making ∆pi p⊥,i p||,i = 0 ror and the ion-cyclotron (IC) instabilities (Hasegawa
0 (where p⊥,i and p||,i are the pressur≡e comp−onents6 of 1969; Gary 1992; Southwood & Kivelson 1993). The
4
ions perpendicular and parallel to B~). Examples of mirror instability consists of non-propagating, strongly
1
compressional modes. Their fastest growing wave
: systems where ion pressure anisotropies are important
v are low-luminosity accretion disks around compact vectors~k are oblique to B~, with the magnetic variations
Xi objects (Sharma et al. 2006), the intracluster medium parallel to B~ being much larger than the perpendic-
(ICM) (Schekochihin et al. 2005; Lyutikov 2007), ular fluctuations, δB~ δB~ . For a bi-Maxwellian
r || ⊥
a and the heliosphere (see, e.g., Hellinger et al. 2006; distribution of ions, and≫assuming cold electrons, the
Matteini et al. 2007; Maruca et al. 2011). In these
threshold condition for mirror instability growth is
systems the ion pressure anisotropy is expected to play
given by T /T 1 > 1/β , where T (T ) is
⊥,i ||,i ⊥,i ⊥,i ||,i
akey roleinthe largescale dynamicsofthe plasma. For −
the ion temperature perpendicular (parallel) to B~, and
instance, Sharma et al. (2006) pointed out that pressure
anisotropiescangiverisetoananisotropicstressthatcan β⊥,i 8πp⊥,i/B2 (Hasegawa 1969). The IC instability,
≡
contribute to angular momentum transport and plasma ontheotherhand,consistsoftransverseelectromagnetic
heatinginlow-luminosityaccretiondisks. Thisphysicsis waves, with the fastest growing ~k being preferentially
notincludedinstandardMHDmodelsofaccretiondisks. parallel to B~ (e.g., Anderson et al. 1991). Whether
the IC or mirror instability sets in first depends on
how fast these instabilities grow for a given set of
1Departamento de F´ısica, Facultad de Ciencias F´ısicas y plasma conditions. Gary et al. (1994) showed that for
Ma2tAemsta´rtoincaosm,yUniDveerpsaidrtamdednetChainled;mTarhioeo.rrieqtuicealml e@Adstfir.oupchhyilsei.ccsl T⊥,i/T||,i − 1 = 0.35/β|0|.,4i2, the growth rate of the IC
Center, University of California, Berkeley, CA 94720; instability is γ =10−4ω , where ω is the cyclotron
IC c,i c,i
[email protected] Center and Department of Physics, frequency of the ions. In addition, for γIC/ωc,i ≪1, the
threshold anisotropy depends very weakly on γ .
University of New Hampshire, Durham, NH 03824; IC
[email protected]
2
In the linear regime, the dominant instability will be which perpendicular heating or parallel cooling drives a
determined by which threshold is reached first as B~ is plasma mirror and/or IC unstable, as in, e.g., the solar
amplified. These estimates for the threshold conditions wind. Hellinger et al. (2008) conducted an analogous
imply that, in the linear regime, the IC instability study to ours using hybrid simulations. In their study,
should dominate for β 1, while the mirror instability expansion of the simulation box decreased the mean
i
should dominate for βi ∼ 1. magnetic field, so that on average p⊥,i < p||,i. This
≫ in turn led to the growth and nonlinear saturation of
Although both the mirror and IC instabilities have been the firehose instability. Also, Travnicek et al. (2007)
extensively studied in the linear regime, a complete conducted a similar hybrid simulation study where the
theoryof their nonlinearevolutionandsaturationis still simulation box is expanded along the directions parallel
lacking. In this paper, we are particularly interested in and perpendicular to B~. With this setup, p becomes
⊥,i
the question of how the mirror and IC modes behave larger than p on average,and both the mirror and IC
||,i
after the initial phase of exponential growth. Indeed, in instabilities can grow and reach saturation. Our work
most astrophysical scenarios where the magnetic field is is complementary to these hybrid simulation studies
amplified, the growth of the field occurs on time scales and focuses on the process in which the condition
muchlongerthanthe growthtime ofthe relevantkinetic p > p is achieved by field amplification due to
⊥,i ||,i
instabilities. Therefore, most of the evolution of the ve- shearing plasmas. We also consider higher values of β
i
locity space instabilities happens in a nonlinear regime, 1 80, relevant to accretion disks, the ICM, and the
≈ −
where the conditions typicallyassumedin linearstudies, near-Earthsolar wind. Our calculations are less directly
such as a homogeneous plasma or a bi-Maxwellian applicabletoheliosphericmeasurementsofmirrormodes
distribution of particle velocities, are not necessarily driven unstable by rapid heating through the Earth’s
satisfied. Moreover, if the field is amplified by order bow shock (e.g., Schwartz et al. 1996). This is because
unity or more, a quasi-linear analysis is not applicable. weassumethatthepressureanisotropyisgeneratedona
One of the key questions we are interested in addressing timescale long compared to the ion cyclotron timescale,
in this case is whether the mirror instability saturates which is probably not true at collisionless shocks.
via pitch-angle scattering that violates magnetic mo-
ment conservation (Sharma et al. 2006) or via a nearly Our paper is organized as follows. In 2 we describe
§
µ -conserving rearrangement of the magnetic field, with the simulation set up, emphasizing the key physical
i
fluctuations δB~ B~ (Schekochihin et al. 2008). We and numerical parameters. 3 shows our results and
shall see that| bo|th∼o|f t|hese saturation processes can in 4 presents our conclusions.§ Throughout the paper
§
fact be important. we compare some of our simulation results with linear
theory predictions appropriate to our PIC simulations.
The nonlinear regime of the IC and mirror insta- These have artificially low ion to electron mass ratios
bilities has been studied both theoretically (e.g., (mi/me 1 10). Thelineartheoryresultsarebasedon
≃ −
Schekochihin et al. 2008; Hellinger et al. 2009) and the linear Vlasov solver developed in Verscharen et al.
numerically (Baumgartel et al. 2003; Travnicek et al. (2013).
2007; Califano et al. 2008; Guo et al. 2009). In most
numerical studies, however, the pressure anisotropy is During the completion of this work, Kunz et al. (2014)
treatedasaninitialcondition,whichdecreasesasthein- presented calculations of firehose and mirror saturation
stabilities grow and saturate. Thus, in these approaches in shearing plasmas with very similar results to those
the nonlinear saturation cannot be followed for a time that we present here.
much longer than that of the initial exponential growth.
Moreover, many of the previous numerical studies have
2. SIMULATIONSETUP
focused on one-dimensional simulations.
We use the electromagnetic, relativistic PIC code
Inthispaperwestudythelongterm,nonlinearevolution TRISTAN-MP (Buneman 1993; Spitkovsky 2005) in
of the mirror and IC instabilities using two-dimensional two dimensions. The simulation box consists of a
particle-in-cell (PIC) simulations in driven systems. In square box in the x y plane, containing plasma with
order to self-consistently explore a time much longer a homogeneous mag−netic field B~ initially pointing
0
than the initial exponential growth, we continuously along the xˆ axis. Since we want to simulate a magnetic
induce the growth of ∆pi by amplifying the mean field field that is being amplified in an incompressible way,
< B~ > during the simulation. We concentrate on we impose a shear motion of the plasma so that the
amplification by incompressible plasma motions. Thus, mean particle velocity is ~v = sxyˆ, where s is a shear
−
our simulations impose a shear velocity on the plasma, parameter with units of frequency and x is the distance
which sustains the growth of the magnetic field and along xˆ. From flux conservation, the y-component
maintains an overall positive pressure anisotropy. This of the mean field evolves as ∂ < B >y /∂t = sB0,
−
growth is intended to mimic magnetic field fluctuations implying a net growth of < B~ > . Although we
in a turbulent plasma and MHD-scale instabilities that present simulations resolvin|g the x y| plane only, we
amplify the magnetic field, like the magnetorotational also tried runs where the x z plan−e was resolved. In
instability (MRI; Balbus et al. 1991) or convective thosecasestheisotropization−efficiencywassubstantially
instabilities driven by anisotropic thermal conduction lower. This is because, if the x z plane is resolved,
in dilute plasmas (Balbus 2001; Quataert 2008). We the growing component of B~ is−perpendicular to the
suspect that our results are also relevant to cases in
plane of the simulation. This way the angle between
Continuously Driven Ion Velocity Space Instabilities 3
the relevant wave vectors, ~k, and B~ cannot be 0. fields (∂/∂x) in the shearing coordinates must include a
This over-constrains the wave vectors that are allowed time dependent term that accounts for the evolution of
toexist,artificiallyreducingtheisotropizationefficiency. the x axis.
The last terms in Equations 1 and 2, which are
Animportantphysicalparameterin oursimulationswill proportional to the coordinate y, will be neglected
be the ratio of the initial ion cyclotron frequency to the for the following reasons. In our simulations the box
shear frequency, ω /s. We refer to this as the magneti- size will typically be a few times the Larmor radius of
c,i
zation. In typical astrophysical environments ω s. the ions, R . Thus, sy/c (s/ω )(v /c) (where
c,i L,i c,i th,i
≫ ∼
Due to computational constraints, however, we will use v is the thermal velocity of the ions), which is
th,i
values of ω /s much smaller than expected in reality, much smaller than unity since we are interested in
c,i
although still satisfying ω /s 1. The dependence the regime s/ω 1 and v /c 1. As a result,
c,i c,i th,i
≫ ≪ ≪
ofourresultsontheratioω /swillbecarefullyassessed. the last term in Equation 1 will be much smaller than
c,i
thetermonthelefthandside,especiallysince B~ E~ .
In standard MHD simulations where shear plasma mo- | |≫| |
tions are imposed, shearing periodic boundary condi- The last term of Equation 2, on the other hand, is
tions would be used along the x direction (see, e.g., not necessarily much smaller than the displacement
Hawley et al. 1995). In that case, the flow velocities at current (left hand side term), mainly because we expect
the x-boundaries of the box are matched using Galilean E~ B~ . However, if the characteristic timescale for
transformations of the fluid velocity. In the case of rel- | | ≪ | |
the mirror and IC modes is close to s−1 (below we will
ativistic PIC simulations this can not be done in a self-
check that this is indeed the case), the last term in
consistent way. The reason is that, under a relativistic
Equation 2 will be much smaller than the first term on
change of reference frame, the current J~ transforms in
the righthandside (the ratiobetweenthese terms scales
a way that is inconsistent with simply transforming the
as (s/ω )2(v /c)2), so we choose to neglect it in
particlevelocities.4 Weinsteadimplementshearingcoor- ∼ c,i th,i
the limit s ω . The fact that the neglected term
dinates,inwhichthegridmoveswiththeshearingveloc- ≪ c,i
can still be comparable to the displacement current
ity ~v = xsyˆ. In this new frame the average plasma ve-
− implies that we may be excluding a relevant contri-
locityvanisheseverywhereintheboxandsimpleperiodic
bution to the charge density in Gauss’ law. However,
boundaryconditionsareallowedbothinthexandyaxes.
the physics of interest in this paper is non-relativistic
In these shearing coordinates, Maxwell’s equations gain
and charge separation does not play a role. As a
additional terms (see Appendix A of Riquelme et al.
result, c B~ + sct∂B~/∂y xˆ 4πJ~ (equivalent to
2012), becoming
∇ × × ≈
c B~ 4πJ~ in the non-shearing coordinate system).
∂B~(~r,t) N∇eg×lectin≈g the terms proportional to y is also required
= c E~(~r,t) sB (~r,t)yˆ+
∂t − ∇× − x by the periodic boundary conditions in the y-direction.
The existence of these terms at all is a consequence
∂E~(~r,t) y∂E~(~r,t)
s ct + xˆ and (1) of the incompatibility of the Galilean invariance of the
(cid:16) ∂y c ∂t (cid:17)× shearing box and the Lorentz invariance of the PIC
calculations. For consistency we will also neglect the
∂E~(~r,t) term sExyˆ in Equation 2, which is also comparable to
=c B~(~r,t) 4πJ~ sEx(~r,t)yˆ the displacement current (in the case where E~ evolves
∂t ∇× − − −
on a timescale close to s−1).5
∂B~(~r,t) y∂B~(~r,t)
s ct + xˆ, (2)
(cid:16) ∂y c ∂t (cid:17)× Our set of simulations are summarized in Table 1. The
physical setup is defined by the magnetization param-
where c is the speed of light, and E~ is the electric field. eter, ω /s (the ratio of the initial cyclotron frequency
c,i
In addition to the modifications to Maxwell’s equations, to the shear rate), the ratio between the initial particle
the forces on the particles also acquire an extra term: and magnetic pressures, β = 8πp /B2, the mass
init init 0
ratio between ions and electrons, m /m , and the initial
dp~ ~u i e
=sp yˆ+q(E~ + B~), (3) Alfv´en velocity of the plasma v , relative to the speed
dt x c × A,0
oflightc, wherev B /√4πρandρ is the massden-
A,0 0
where p~, ~u, and q are the particle’s momentum, velocity, sity. The numerical p≡arameters of our runs are defined
and charge, respectively. The third term on the right by the spatial resolution (c/ω /∆ ), box size relative
p,e x
hand side of Equation 1 and the fourth term on the to the ion Larmor radius (L/R ), and number of par-
L,i
right hand side of Equation 2 are proportional to time ticles per cell (N ), where ω is the plasmafrequency
ppc p,e
and arise from the motion of the shearing coordinate of the electrons and ∆x is the spacial separation of the
grid points with respect to those of the non-shearing gridpoints. The simulations shownin this paper arede-
coordinates. Indeed, as time goes on, the x axis of the scribed in Table 1. These simulations are a subset of
shearingcoordinatesisgraduallytilted(asseenfromthe a higher number of runs, where we used different com-
non-shearing frame); therefore the x derivatives of the
−
5 The modification to Faraday’s equation that depends on Bx
4 Relativisticvelocitytransformationsimplychangesinthevol- can be integrated using simple time and space interpolations of
ume occupied by the particles, leading to a non-conservation of Bx. Thisway,afterthismodificationisimplemented,thenumeri-
chargedensitythatcannotbetriviallyimplementedinaPICsim- calalgorithmusedbyTRISTAN-MPcontinuestobesecondorder
ulation(seeRiquelmeetal.2012). accurateintimeandspace.
4
TABLE 1
Physicalandnumericalparametersforthe runs
Runs mi/me βinit βinit,i (=βinit,e) ωc,i/s vA,0/c c/ωp,e/∆x Nppc L/RL,i
beta6mag93 1 6 - 93 0.05 14 30 12
beta6mag670a 1 6 - 670 0.15 14 10 18
beta6mag670b 1 6 - 670 0.15 14 10 24
beta20mag93a 1 20 - 93 0.05 10 60 35
beta20mag93b 1 20 - 93 0.05 10 30 35
beta20mag670a 1 20 - 670 0.05 7 30 36
beta20mag670b 1 20 - 670 0.05 7 30 48
beta20mag2000 1 20 - 2000 0.05 7 50 55
beta80mag670 1 80 - 670 0.05 5 30 24
beta80mag1340 1 80 - 1340 0.05 5 30 24
betai20magi240mass10 10 - 20 240 0.05 5 20 25
Note. — A summary of the physical and numerical parameters of the simulations discussed in the paper.
Thesearethemassratiomi/me,βinit (=βinit,i+βinit,e)ifmi/me=1,βinit,i (=βinit,e)ifmi/me>1,the
magnetization ωc,i/s,theinitialAlfv´envelocity, vA,0,theskindepth c/ωp,e/∆x (where∆x isthegridpoints
separation), the number of particles per cell Nppc (including ions and electrons), and the box size in units of
the typical ionLarmor radiusL/RL,i (RL,i =vth,i/ωc,i, wherevt2h,i =3pi/ρis the rmsionvelocity and ρis
themass densityof the ions). We confirmed numericalconvergence byexploringtheresolutioninc/ωc,e/∆x,
Nppc,andL/RL,i forallparametercombinations ofβinit (orβinit,i),ωc,i/s,andmi/me.
Fig.1.—ThethreecomponentsofδB~ andplasmadensityfluctuationsδρattwodifferenttimes: t·s=1(upperrow)andt·s=2(lower
row), for run beta6mag670b (βinit = 6, ωc,i/s =670, mi/me = 1). Fields and density are normalized by B0 and the initial density ρ0,
respectively. Arrowsdenotethemeanmagneticfielddirectiononthesimulationplane. ICandmirrorinstabilitiescontributecomparablyto
thefluctuationsatt·s=1,whilemirrordominatesatt·s=2. Densityfluctuationsbestcorrelatewiththemirrormodes. Thedominance
of mirror fluctuations suggests that mirror modes are more robust than the IC modes in the saturated regime. This may be due to the
particleenergyspectrumdepartingfrombi-Maxwellian(seeFigure5),suppressingthegrowthoftheICmodes.
binations of numerical and physical parameters, which and β = β + β to quantify the pressure of the
|| e,|| i,||
confirmed numerical convergence. plasma perpendicular and parallel to B~. Simulations
with m /m = 1 of course do not allow us to study
i e
the physics of electron isotropization. However, in 3.4
3. RESULTS we use an example of our mi/me > 1 runs to sh§ow
that m /m = 1 calculations reproduce the ion-related
We want to explore the regime β ,β 1 80. In i e
⊥ ||
≈ − phenomenafairlywell, withβ andβ playingthe same
order to do so, most of our runs use the mass ratio ⊥ ||
mi/me = 1. Since for mi/me = 1 both species will be role of βi,⊥ and βi,|| in the mi/me > 1 runs. We defer
essentially indistinguishable, we initially give ions and a detailed study of electron scale pressure isotropization
electronsMaxwellianenergydistributions with the same to future work.
temperature, and use the parameters β = β +β
⊥ e,⊥ i,⊥
Continuously Driven Ion Velocity Space Instabilities 5
Fig.3.—Evolutionoftheionmagneticmomentfortheβinit=6
runswithωc,i/s=93and670(beta6mag93 andbeta6mag670a in
Table1). Weshowtheevolutionofthetrueaveragemagneticmo-
mentµ≡<v⊥2/B>p (blackline)andaneffectiveglobalmagnetic
moment, µeff ≡< v⊥2 >p /| < B~ > | (red line), where the sub-
script p denotes an average over all particles. The ion magnetic
momentµisconserveduntilts∼1,atwhichpointitdecreaseson
the same timescale as the mean magnetic field (see Fig. 4). The
closesimilaritybetweenµandµeff showsthattherearenotstrong
correlationsbetweentheparticlevelocityv⊥andmagneticfieldB;
suchcorrelationsaremoreprominentathigherβinit wheretheIC
instabilityissubdominant(seeFig. 9).
Fig. 2.— Field fluctuations and pressure anisotropy for our
βinit = 6 runs. Panels a and b show the evolution of δBj2/B02
(≡<(Bj−<Bj >)2 >/B02; solid) and Bj2/B02 (≡<Bj >2 /B02;
dotted)forrunsbeta6mag93andbeta6mag670a,respectively(j=
x,y,andz correspondtoblack, red,andgreen, respectively). For
the sameruns,panels cand dshow ∆p/p|| (≡<(p⊥−p||)/p|| >;
black-solid), compared with the linear mirror and IC thresholds
(∆p/p and∆p/p ,inredandgreen,respectively)forpair
||,MI ||,IC
plasma and growth rates γ =0.05ωc,i and γ =0.007ωc,i, respec-
tively. Panels e and f show δB2/B2 and δB2/B2, where the
|| ⊥
subscripts||and⊥denotethecomponentsparallelandperpendic-
ularto<B~ >,respectively. Aftertheplasmapressureanisotropy Fig.4.—Therateatwhichtheionmagneticmomentchangesin
exceeds the mirror and IC thresholds, there is an initial phase time,|dlnµ/dt|,fordifferentβinit,comparedtothegrowthrateof
of exponential growth followed by a secular phase. The pressure themeanmagneticfielddlnh|B~|i/dt(blackline;thisquantityisthe
anisotropyinitiallygrowsbutsaturates atthelinearthresholdfor
themirror/ICinstabilities. sameforthedifferentβinitruns). |dlnµ/dt|isshownforoursimu-
lationswithβinit=6,20,and80(beta6mag670a, beta20mag670a,
andbeta80mag1340, respectively). Thevertical-dotted linesmark
We divide our runs into cases with three different initial the beginning of the saturated state for each run based on when
betas: β =6, 20,and80. Our analysisfocuses onthe thegrowthofthefluctuationsbecomesroughlysecular. Inthesat-
init
nonlinear structure of the IC/mirror generated fluctua- uratedstate,themagneticmomentchangesonatimescalecompa-
rable to that of the mean magnetic field (to within ∼25−50%).
tions,theevolutionofthepressureanisotropyp /p 1,
⊥ || Atearliertimes,however,themagneticmomentisreasonablywell
−
and the conservation of the ion magnetic moment. conserved.
3.1. Case β =6 for the component x, y, or z of δB~) at two different
init
times: t s=1(upperpanels)andt s=2(lowerpanels).
As discussed above, the expectation is that the IC · ·
instabilityshouldplayanimportantroleinthe β =6
init We can see that at the earlier time (t s = 1) the three
case, and it should become significantly less important ·
components of δB~ have about the same amplitude (see
in the β = 20 and 80 cases. The contributions
init Figures 1a, 1b, and 1c), but appear to be dominated by
from the IC and mirror modes can be seen from
different mode orientations, indicating the simultaneous
Figure 1. This figure shows 2D images of the field
presence of IC and mirror modes. The mirror modes
fluctuations δB /B ( (B < B >)/B ; <>
stands for “volujme0ave≡rage”)j−and the jplasma 0density are expected to have δB~ mainly in the plane of the
fluctuations δρ/ρ for run beta6mag670c (“j” stands simulation (the plane of~k and <B~ >; Pokhotelov et al.
0
6
2004). On the other hand, the IC modes have the three
components of δB~ of comparable magnitude, so their
presence can be most clearly revealed by δB (Figure
z
1c). Indeed, while δB and δB are dominated by
x y
a combination of oblique (mirror) and quasi-parallel
(IC) modes, δB appears to be dominated by only
z
quasi-parallel modes. Thus, initially the mirror and
IC instabilities contribute comparably to δB~. At the
later time (t s = 2), however, the fluctuations become
·
dominated by an oblique wave vector, with the IC
instability playing a subdominant role. The density
fluctuations at all times are δρ/ρ 1, and seem to
0
≪
correlate primarily with the mirror modes.
The relative contribution of the mirror and IC in-
stabilities can also be seen from Figure 2, which
shows the time evolution of different volume-averaged
quantities for runs with ω /s = 93 (first column;
c,i
run beta6mag93a) and ω /s = 670 (second column;
c,i
run beta6mag670a). In Figures 2a and 2b we plot
δB2/B2 for the two magnetizations. In both cases
j 0
there is an initially exponential growth of δB~ until
| |
< δB2 > /B2 0.03. Both instabilities have about
j 0 ∼
the same growth rates, which are γ γ 10s.
IC MI Fig.5.— The energy spectra of particles for the high magneti-
≈ ∼
Note that the dominant growth rate of the mirror/IC zation βinit = 6 run beta6mag670a at two different times. Top:
modes is nearly the same in the two simulations with Energyspectraasafunctionoftheenergiesperpendicularandpar-
allelto the magnetic field, p2/2m (black-dotted line) and p2/2m
different values of ωc,i. The growth rate is thus set by ⊥ ||
the background shear rate s instead of ωc,i. This can (black-solidline). Bottom: Energyspectraasafunctionofvk. Bi-
be understood by noting that the smaller growth rate Maxwellianfits are shown in red. There is a clear deviation from
the bi-Maxwellian at the highest energies, which grows in time.
modes have a smaller anisotropy threshold. Therefore,
There is, however, little deviation from a Maxwellian at v ∼ 0,
k
as B and ∆p grow, modes with smaller growth rate will in contrast to some of the 1D runs discussed in the Appendix.
begin to grow first. This implies that the dominant Comparison with Figure 10 suggests that the high energy devia-
modes will be those that can reach a significant ampli- tion from Maxwellian is due to the IC instability, since it is less
tudetostopthegrowthof∆ponthesheartimescales−1. prominentathigherβ whentheICinstabilityissub-dominant.
the isotropizationoccurs.
After the exponential phase, the growth becomes domi-
natedbythemirrormodesandissecular(aspredictedby
One consequence of the difference in δB~ is that the
Schekochihin et al.2008). Thistransitionfromexponen- | |
lower magnetization runs take longer to reach the point
tial to secular is clearly seen at ts 1 for ω /s = 670
in Figure 2. Despite the similar exp∼onential gcr,iowth rate when ∆p/p|| (≡< (p⊥ −p||)/p|| >) saturates. Thus in
the ω /s=93runthe anisotropy∆p/p hasmoretime
for the mirror and IC instabilities in the linear regime, c,i ||
to grow, reaching larger values before saturating. This
the mirror modes are more robust and dominate in the
can be seen in the evolution of ∆p/p , shown in black
nonlinearregime. Thisdominanceofthemirrormodesis ||
lines in Figures 2c and 2d for runs beta6mag93a and
remarkable because it occurs even for β 1, where the
i
∼ beta6mag670a, respectively. For both magnetizations,
linear analysis predicts an important contributionof the
there is an initial overshoot in ∆p/p , whose amplitude
IC instability (in the nonlinear regime the mean mag- ||
netic energy has been amplified by about one order of is larger for the ωc,i/s = 93 run. After the overshoot,
magnitude, so that βi 1). We suggest below (Fig. 5) ∆p/p|| evolves in a similar way for both magnetizations.
∼
that departure from the bi-Maxwellian energy distribu-
tion assumed in the linear instability analysis is likely Sinceinthe nonlinearregimeδB~ isdominatedby mirror
playing an important role in this sub-dominance of the modes, after the saturation one would expect the pres-
IC instability. sure anisotropies to behave according to the marginal
stability conditions of the mirror instability. Since in
OnedifferenceintheevolutionofδB~ forthetwodifferent our simulations both the IC and mirror modes grow at
magnetizations is that the end of the exponential phase γ 10s, we calculated the IC and mirror thresholds
∼
occurs at somewhat higher amplitude in the ω /s=93 for equivalent growth rates in the pair plasma case with
c,i
case. This difference implies that, for smaller values of magnetizations ωc,i/s = 93 and ωc,i/s = 670. Thus,
ω /s, the mirror and IC fluctuations are less efficient Figures 2c and 2d show that the threshold for mirror
c,i
in suppressing the pressure anisotropy, requiring larger (red) and IC instability (green) in the cases γ =0.05ωc,i
values of δB~ . At the end of this section, we will and γ = 0.007ωc,i, respectively. We see that the mirror
| | threshold coincides fairly well with the saturated ∆p/p
explain in further detail this dependence of δB~ on ||
| | for both magnetizations. These results imply that the
magnetization by focusing on the mechanism by which
linear marginal stability condition for the mirror insta-
Continuously Driven Ion Velocity Space Instabilities 7
µ corresponds to the actual average of the magnetic
moment. Thus if a nearly µ conserving process dom-
−
inates the mirror saturation, µ should remain fairly
constant. µ , on the other hand, can be thought as
eff
an effective global magnetic moment that ignores the
fluctuations in B~ or the correlation between v and B
⊥
due to particles collecting in mirrors. µ necessarily
eff
decreases since the mirror and IC instabilities suppress
the growth that < v2 > would have if the particles
⊥ p
were only affected by the mean field < B~ >. Thus, if
µ is nearly conserved, µ and µ should evolve quite
eff
differently, with µ > µ . On the other hand, if there
eff
is significant pitch angle scattering µ µ and both
eff
∼
shoulddecreaseonthesametimescalethat B~ increases.
h i
Figures 3a and 3b show µ and µ for our β = 6
eff init
runs with magnetizations ω /s = 93 and ω /s = 670,
c,i c,i
respectively. We see that for both magnetizations the
differencebetweenµ andµ isverysmall(only a few
eff
∼
% difference). This implies that there is little spatial
Fig.6.— The time derivative of the total particle thermal correlation between v and B. We shall see below that
energy(black)alongwiththevolumeaveragedheatingratebythe ⊥
µ and µ differ somewhat more at higher β where the
anisotropic stress −s∆pˆbxˆby (red) for βinit =6 and ωc,i/s=670 mirror ineffstability dominates over the IC instability.
(runbeta6mag670a). Thegoodagreementbetweenthetworesults
implies that the anisotropic stress coupling to the background
shearistheprimarymechanismforparticleheating. Figure4comparestherateofchangeoftheionmagnetic
moment with that of the mean magnetic field for runs
with β =6, 20,and 80. This comparisonis important
init
bilityisfairlyaccurateindetermining∆p/p inthenon-
|| since it is the evolution of the mean field that is driving
linear regime. In realistic astrophysical setups, where
the velocity space instabilities in our calculations. If the
ω can be many orders of magnitude larger than s, the
c,i totalthermalenergyvariesonatimescalelongcompared
expectation is that ∆p/p will follow the γ 0 mirror
|| to the mean magnetic field, maintaining marginal sta-
→
threshold, givenby p⊥/p|| 1=1/β⊥. This differs from bility to the mirror instability via pitch-angle scattering
−
theresultsshowninFigure2,whichareonlyappropriate implies dlnµ/dt dlnB/dt[1 (β−1)]. The term
for pair plasmas and modest magnetization. (β−1)|is exactl|y≃β−|1 for th|e c−anOonical high magne-
O
tization mirror instability threshold but is somewhat
Asdiscussedintheintroduction,itisnotknownwhether
different in our pair plasma, modest magnetization
themirrorinstabilitysaturatesviapitch-anglescattering
simulations (hence the use of ). This expression for
(see, for instance, Sharma et al. 2006) or whether it O
dlnµ/dt implies that pitch-angle scattering should lead
cancels the appearance of a pressure anisotropy by
to the magnetic moment varying at a rate slightly less
substantially modifying the structure of the magnetic
than that of the mean magnetic field. This is consistent
field (Schekochihin et al. 2008). In the latter case δB~ with our numerical results in Figure 4 in the saturated
should grow secularly until δB~ B, with the breaking state where the magnetic field grows secularly in time.
| |≃
of µ invariance not necessary for the regulation of p .
⊥
Figu−res 2e and 2f show the evolution of δB2/B2 and At earlier times, when the magnetic field fluctuations
⊥
δB2/B2 for ω /s=93 and 670, respectively, where B~ due to the mirror and IC instability are smaller, µ is
|| c,i ⊥ approximatelyconservedeventhoughthe mirrorand IC
andB~|| refertothemagneticfieldcomponentperpendic- are present and grow exponentially. In particular, there
ular and parallel to the mean field < B~ >. We see that is a temporal lag of 0.3s−1 between the onset of the
∼
i) δB2 dominates in the nonlinear regime, as expected exponential growth of the mirror and IC instabilities
|| and the onset of pitch angle scattering that decreases µ.
from the larger amplitude of the oblique mirror modes,
Our interpretation of this result is that the mirror/IC
and ii) the saturation amplitude is δB~ 2/B2 0.04,
fluctuations must reach a sufficient amplitude, roughly
| | ∼
which favors the scenario where the non-conservation
δB 0.1 0.3 B , in order for pitch angle scattering to
of µ is the key mechanism for the isotropization of the ∼ − h i
be effective. For the mirror instability, this violation of
plasma pressure, at least if the background field has
µ conservation by finite amplitude fluctuations is likely
been amplified significantly.
due to the stochasticity of particle orbits that sets in
for large amplitude (albeit low frequency) fluctuations
To quantify the time variation of the average ion mag-
(e.g., Chen et al. 2001; Johnson et al. 2001). For the
netic moment we define
ion cyclotron instability, the need for finite amplitudes
v2 <v2 > before significant scattering sets in is a consequence
µ=(cid:28)B⊥(cid:29) and µeff = <B⊥~ >p , (4) of the well-known result that the scattering rate for
p | | cyclotron-frequency fluctuations is ωc,i(δB/B)2.
∼
where <> stands for average over all the particles.
p
8
Fig.7.—Magneticfieldfluctuations andpressureanisotropyforourβinit=20runswithωc,i/s=93,670,&2000 (runsbeta20mag93a,
beta20mag670a,andbeta20mag2000). Theresultsforωc,i/s=670and2000arequantitativelysimilar,indicatingthatourresultsaccurately
describe the high cyclotron frequency limitrelevant to astrophysical and heliospheric systems. Panels a-c show the evolution of δBj2/B02
(≡<(Bj−<Bj >)2>/B02;solid)andBj2/B02 (≡<Bj >2/B02;dotted),respectively(j=x,y,andzcorrespondtoblack,red,andgreen,
respectively). For the same runs, panels d-f show ∆p/p|| (≡< (p⊥−p||)/p|| >; black-solid), compared with the linear mirror threshold
forapairplasmaandgrowthrates γ=0.05ωc,i (panel d)andγ=0.007ωc,i (panels eandf). Panels g-ishowδB|2|/B2 andδB⊥2/B2. In
contrast with the βinit =6case (shown inFigure2), here δB|2| ≫δB⊥2, which isconsistent withthe dominance of the mirrorinstability.
Alsonotethatthenonlinearsaturationofthemirrorfluctuations occurswhen|δB~|∼0.3B.
The breaking of local µ-invariance for this calculation bi-Maxwellian for our β = 6 run, mainly due to a
init
(β , β < β = 6) might be affected by the subdomi- significant excess of particles at the highest energies.
⊥ || init
nant (but still significant) contribution of the IC modes In future work, we will study whether or not the
to the magnetic fluctuations. However, Figure 4 shows spectral variation in Figure 5 can completely account
thatthemagneticmomentvariesatthesamerateasthe for the decrease in the amplitude of the IC modes in
meanmagneticfieldforalloftheβ wehavesimulated. the saturated state. Figure 5 also shows that the ion
init
We will discuss the higher β results in more detail in velocity distribution does not deviate significantly from
init
the next two subsections. abi-Maxwellianatlowvk,ashasbeenfoundinsome1D
studies of the saturation of the mirror instability (e.g.,
An interesting question is whether the non-linear evolu- Southwood & Kivelson 1993; Califano et al. 2008). We
tion of the velocity space instabilities makes the energy discuss this in more detail in the Appendix.
distribution of the particles substantially different from
a bi-Maxwellian spectrum. This question is particularly Finally,Figure6showstheevolutionofthetotalheating
relevant in terms of understanding the suppression of rateoftheparticlesasafunctionoftime(blackline). For
the IC modes in the nonlinear regime as β goes from comparison, in the red line we plot the particle heating
⊥
6 to 1. Indeed, it has been shown(Isenberg 2012; rate by the anisotropic stress, s∆pˆb ˆb (Sharma et al.
Isenbe∼rg et al. 2013) that a significant departure from a 2006). We see a reasonable ag−reemenxt byetween the two
bi-Maxwellian distribution can increase the anisotropy results, implying that the anisotropic stress plays the
threshold for the growth of the IC modes and that dominant role in particle heating. Other mechanisms,
a change of this kind is expected given the resonant likewave-particleinteraction,mustplayasecondaryrole
character of the IC instability. Figure 5 shows that in the total particle energization. This is true for all of
the particle spectrum does differ somewhat from a the β we have simulated.
init
Continuously Driven Ion Velocity Space Instabilities 9
Fig.9.—Panelsaandbshowtheevolutionoftheionmagnetic
moment µ for βinit = 20 runs with ωc,i/s = 93 (beta20mag93a;
panel a) and ωc,i/s = 670 and 2000 (beta20mag670a and
beta20mag2000; panel b). Panels c and d compare the corre-
Fig.8.—SpatialdistributionofmagneticfieldfluctuationsB2/< sponding rate of change of µ with that of the mean magnetic
B2 > at t·s ≈ 1 for βinit = 20 and magnetization ωc,i/s = fieldforrunsbeta20mag93a andbeta20mag670a, respectively. All
9b3eta(p2a0nmealga6;70rbu)n.bPeatnae2l0smcaagn9d3bd)shaonwdtωhce,is/asm=eq6u7a0nt(iptaiensealtbt;·sru≈n p≡l<otsv⊥s2h/oBw>thpe(ebvloalcuktiolinneo)fatnhde tarnueeffaevcetriavgeegmloabganlemticagmnoetmicenmtoµ-
2. The black arrows show the direction of < B~ >. The higher ment, µeff ≡< v⊥2 >p /| < B~ > | (red line), where the subscript
magnetization run produces mirror modes whose wavevectors are pdenotes anaverageoverallparticles. Themagneticmomentbe-
more perpendicular to < B~ > and with longer wavelengths. A ginstodecreasewhenthefluctuationsreachδB∼0.1hBi(seeFig.
migrationtolargerwavelengths intimecanalsobeseen. 7). Thedifference between µandµeff ismoresignificantthan for
βinit=6(cfFig. 3). Thisisaconsequenceofparticlesbunchingin
3.2. Case βinit =20 mlesirsr,otrhse,wiohnicmhalegandestictomaocmorernetlacthioanngbeestwateeanbovu⊥tatnhde Bsa.mNeornaettehaes-
We now consider a somewhat weaker magnetic field thebackgroundmagneticfieldaftersaturationatts≃0.6,indicat-
case with β = 20. In this case, the field fluctuations ing that pitch angle scattering regulates the nonlinear saturation
init onsufficientlylongtimescales.
are dominated by the mirror instability at all times.
This can be seen in Figures 7a-c which show the
evolution of the three components of δB~2 for β = 20 the mirror’s dominant wave vector, ~k, with respect to
runs with magnetizations ω /s=93, 670, anindit2000, < B~ >: a large δB~2/δB~2 ratio implies that ~k and
c,i || ⊥
respectively (runs beta20mag93a, beta20mag670a and < B~ > are nearly perpendicular. The results of Figures
beta20mag2000a in Table 1). In all three cases δBz2 7g-i imply that for larger magnetization, ~k and < B~ >
is small compared to δB~2, showing the subdominant are more perpendicular, which is consistent with linear
role of the IC instability. Figures 7a-c also show the calculations (Pokhotelov et al. 2004).
characteristic transition between exponential growth
(with γ 10s) and secular growth (with δB~2/B2 0.1) The orientation of the modes in the nonlinear stage can
≈ ∼
for the mirror modes. alsobe seendirectlyinFigure8,whichshowsthe spatial
distribution of B~2/ B2 at two different times for the
As in the βinit = 6 case, the transition between the runs with ωc,i/s = 9h3 aind ωc,i/s = 670. Figures 8a and
exponential and secular regimes occurs at smaller 8b correspondto the time t s=1, while Figures 8c and
amplitudes and in a smoother way in the more strongly 8d correspond to t s = 2. ·For each magnetization, the
magnetized runs. Figures 7d-f also show that, as in mirror fluctuations·are dominated by two modes that
the case of βinit = 6, after the overshoot the pressure are symmetric with respect to the magnetic field, and
anisotropy follows the marginal stability condition for that are more oblique at higher magnetization. Note
the mirror modes reasonably well, particularly at higher also that the wavelength of the dominant modes (in
magnetization. units of the Larmor radius of the particles, R ) is
L,i
larger for larger magnetization and that the modes tend
Figures 7g-i show δB~2/B2 and δB~2/B2 as a function to grow in wavelength as time goes on.
|| ⊥
of time. While δB~2/B2 is about the same for the three
|| Figure9showstheevolutionoftheionmagneticmoment
magnetizations, δB~2/B2 is significantly smaller in the µandtherateofchangeofµrelativetothatofthemean
⊥
higher magnetization runs. The relative magnitude magneticfieldforourβ =20calculations. Acompar-
init
of δB~ and δB~ is a measure of the orientation of ison of Figures 9 & 7 shows that the magnetic moment
|| ⊥
10
gardless of the magnetization. Since these simulations
are almost purely mirror dominated, this suggests that
themirrormodesproviderathermomentum-independent
pitch-angle scattering to the particles. As a result, it
is likely that IC modes produce the deviations from a
bi-Maxwellian distribution seen in Figure 5. Figure 10
also shows that the parallel velocity spectrum dN/dv
k
remains nearly Maxwellianduring most of the saturated
regime (with only a transient 1% decrease for v 0
||
∼ →
especially at the end of the exponential growth regime).
The lack of significant flattening of dN/dv (which con-
k
trasts with the prediction of Califano et al. 2008) is an-
alyzed in further detail in the Appendix.
3.3. Case β =80
init
For β = 80 we focus on simulations with ω /s =
init c,i
670 and 1340, respectively (runs beta80mag670 and
beta80mag1340 in Table 1). The results are essentially
the same as for the β =20 case so we do not discuss
init
them in detail. Figure 11 shows that the mirror insta-
bility dominates the fluctuations in B~. The evolution of
the pressure anisotropy ∆p/p contains the initial over-
||
shoot followed by a rather flat behavior. The amplitude
oftheovershootiscontrolledbythemagnetization,while
Fig. 10.—Energyspectraofparticlesforourhighmagnetization the subsequent evolution is fairly well described by the
βinit = 20 run beta20mag670a at two different times during the threshold condition for mirror modes.
nonlinearregime. Top: Energyspectraasafunctionoftheenergies
perpendicular and parallel to the magnetic field, p2/2m (black-
⊥
dottedline)andp2/2m(black-solidline). Bottom: Energyspectra As in the βinit =20 case, the saturation δB/B is 0.3.
|| ∼
Figure 12 shows the evolution of the ion magnetic
asafunctionofv . Inallcases,abi-Maxwellianenergyspectrum
k
(red)providesagoodapproximationtothenumericalsimulations. moment for the βinit = 80 simulations. The evolutions
of µ and µ are the same for the two magnetizations
eff
is reasonably conserved until the magnetic fluctuation presented (ω /s = 670 and ω /s = 1340), and show
c,i c,i
amplitude is close to (though somewhat less than) its the same behavior seen in the β = 20 case, which is
init
saturated value. As in our βinit = 6 calculations, this discussed in detail in the previous subsection.
implies that there is a time lag of 0.25 0.3s−1
∼ −
betweenthe onsetofthe mirrorinstabilityandthe onset
of significant pitch-angle scattering (and its associated 3.4. Comparison with m /m >1 simulations
i e
decrease in µ). Note, however, that for β = 20 the
init In this section we show that our use of m /m = 1
saturated amplitude of the fluctuations in the secular i e
in the previous sections does not have a significant
regime, δB~ 0.3B, is larger than in the βinit = 6 effect on the evolution of the ion pressure anisotropy
| | ∼
case, for which δB~ 0.1B. This implies that when under the influence of the IC and mirror instabilities.
only the mirror|ins|ta∼bility plays a significant role, a In Figure 13 we present the results for a simulation
somewhat larger fluctuation amplitude is necessary for with m /m = 10, β = β = 20, and ω /s = 240
i e i e c,i
efficient pitch-angle scattering. We will see below that (runbetai20magi240mass10c),whichis analogousto the
δB~ 0.3B continues even for βinit =80; the saturated previous mi/me = 1 simulations that use βi = βe = 10.
a|mp|l∼itudeofthemirrormodesisthusfairlyindependent A comparison of Figure 13a with Figures 7a, 7b, and
of the beta of the plasma as long as the IC instability is 7c shows that for both mass ratios the evolution of δB~
not significant (β & 6). is dominated by the mirror instability, with the initial
exponential regime being followed by secular growth.
Figure 9 also shows that for β = 20, the true ion Figure13b showsthe evolutionof δB2/B2 andδB2/B2.
init ⊥ ||
magnetic moment µ decreases somewhat more slowly As in the m /m = 1 case (see Figures 7g, 7h, and 7i),
i e
than the effective global magnetic moment µeff (by δB2 δB2 in the secular stage, implying that mirror
25%). This was not seen in our β =6 calculations || ≫ ⊥
(∼Figs. 3). The modest difference ibneittween µ and µeff modes with k⊥ ≫ k|| are dominant. The maximum
amplitude of δB /B 0.3 is also in good agreement
indicates that the nonlinear saturation of the mirror ||
∼
instability involves correlations between v⊥ and B, i.e., with the mi/me =1 simulations.
particles bunching in mirrors. This effect is independent
of magnetization, as can be seen from the solid and The evolution of the ion pressure anisotropy ∆pi/p||,i is
dotted lines in Figure 9b (corresponding to ω /s= 670 shown in Figure 13c. As in the m /m = 1 runs (see
c,i i e
and 2000, respectively). Figures 7d f), ∆p /p follows the linear threshold of
i ||,i
−
themirrorinstabilityfairlyclosely. ThisisseeninFigure
Figure 10 shows that the particle energy spectrum stays 13c,whichshowsthelinearmirrorthresholdexpectedfor
fairly close to bi-Maxwellian in our β = 20 runs, re- modes with γ/ω = 0.05, β = β = 20, and m /m =
init c,i i e i e