Cosmic Ray propagation in sub-Alfvenic magnetohydrodynamic turbulence

0.68 MB·

by  Romain Cohet

Astronomy&Astrophysicsmanuscriptno.AA-27376 (cid:13)cESO2016
March25,2016
Cosmic ray propagation in sub-Alfvénic magnetohydrodynamic
turbulence
R.Cohet1 andA.Marcowith1
LaboratoireUniversetparticulesdeMontpellier,UniversitédeMontpellier,CNRS/IN2P3,CC72,PlaceEugèneBataillon,
34095MontpellierCedex5,France
e-mail:[email protected]
Received-;accepted-
6
1 ABSTRACT
0
2 Context.Thepropagationofcosmicraysorenergeticchargedparticlesinmagnetizedturbulenceisacomplexproblemwhichinvolves
non-linearwave-particleinteractionsandchaoticmagneticfieldlinestransport.Thisproblemhasbeenaddresseduntilrecentlyusing
r
a eitheranalyticalcalculationsorsimulationsusingprescribedturbulencemodels.Withtheadventofsupercomputersitisnowpossible
M toinvestigateenergeticchargedparticlepropagationusingdirectcomputationofelectromagneticfields.Thisisinparticularthecase
forhigh-energyparticlespropagationinmagnetohydrodynamicturbulence.
Aims. This work has the main objective to provide a detailed investigation of cosmic ray propagation in magnetohydrodynamic
4
turbulent fields generated by forcing the fluid velocity field at large scales. It provides a derivation of the particle mean free path
2
dependencesintermsoftheturbulenceleveldescribedbytheAlfvénicMachnumberandintermsoftheparticlerigidity.

Methods. We use an upgrade version of the magnetohydrodynamic code RAMSES which includes a forcing module and a kinetic
]
E moduleandsolvetheLorentzequationforeachparticle.Thesimulationsareperformedusinga3dimensionperiodicalboxinthe
H test-particleandmagnetostaticlimits.TheforcingmoduleisimplementedusinganOrnstein-Uhlenbeckprocess.Anensembleaverage
overalargenumberofparticletrajectoriesisappliedtoreconstructtheparticlemeanfreepaths.
h. Results.WederivethecosmicraymeanfreepathsintermsoftheAlfvénicMachnumbersandparticlereducedrigiditiesindifferent
p turbulenceforcinggeometries.Thereducedparticlerigidityisρ=rL/LwhererListheparticleLarmorradiusandListhesimulation
- box length related to the turbulence coherence or injection scale Linj by L ∼ 5Linj. We have investigated with a special attention
o compressibleandsolenoidalforcinggeometries.
r Conclusions.Wefindthatcompressibleforcingsolutionsarecompatiblewiththequasi-lineartheoryormoreadvancednon-linear
t
s theorieswhichpredictarigiditydependenceasρ1/2orρ1/3.SolenoidalforcingsolutionsatleastatlowormoderateAlfvénicnumbers
a arenotcompatiblewiththeabovetheoreticalexpectationsandrequiremorerefinedargumentstobeinterpreted.Itappearsespecially
[ forAlfvénicMachnumbersclosetoonethatthewanderingoffieldlinescontrolsperpendicularmeanfreepathsolutionswhatever

theforcinggeometry.
2
v Key words. Physicaldataandprocesses:Magnetohydrodynamics(MHD)-Turbulence–Methods:numerical–Interstellarmedium:
1 CosmicRays
7
9
4
1. Introduction field) diffusion coefficients at 90 degrees and hence produces a
0
. pathological behavior in the calculation of the particle parallel
1
Cosmic rays (CRs) propagate in our Galaxy through the in- mean free path (mfp). Several analytical attempts have been
0
teraction with electromagnetic perturbations usually described proposed to deal with this issue (see Shalchi (2009) and refer-
6
in the magnetohydrodynamic (MHD) approximation (Schlick- ences therein). In particular, the 90 degrees scattering problem
1
:eiser 2002), i.e. long-wavelength perturbations with scales hasbeentreatedbythemeanofthebroadeningoftheresonance
vcomparable to the particle’s Larmor radius r = vγ/Ω , where (Völk1973).However,thewaytheresonanceisbroadenedand
L c
XiΩc =ZeB/mcisthecyclotronpulsationofaparticleofvelocity the non-linear solutions which have to be retained in diffusion
v, mass m, and charge Z in a magnetic field of strength B. The coefficient calculations, are model dependent. Needless to
r
amodeling of particle transport in MHD turbulence is complex. say that all these studies are performed in the test-particle
It requires approximations to obtain the solutions that describe approximation(seePtuskinetal.(2006)).Itisknownthatifthe
theparticletrajectories.Transportstudieswereoriginallybased density in the relativistic particles is high enough, the particles
on the quasi-linear theory (QLT), (Jokipii 1966; Schlickeiser can generate their own waves and induce a modification of the
2002)wheretheunperturbed(gyromotion)particletrajectoryis turbulence (Farmer & Goldreich 2004; Yan & Lazarian 2008,
retained to derive the electromagnetic correlation tensor, itself 2011).
usedtoderivetherandomLorentzforceexertedontheparticles. AnotherdifficultyinCRtransportstudiesliesinthedescription
The QLT is applicable over restricted timescale ranges, i.e. of the random Lorentz forces. In the incompressible MHD
intermediatebetweenpitch-angleparticlescatteringandparticle limit the turbulence can be described by an anisotropic model
distribution isotropization times. The main drawback of the (Goldreich & Sridhar 1995; Maron & Goldreich 2001). In this
QLT is that it leads to infinite cosine pitch-angle (the angle model, owing to the Alfvén wave packet dynamics, the turbu-
between the particle velocity and the background magnetic lent perturbations are elongated along the magnetic field and
Articlenumber,page1of15
A&Aproofs:manuscriptno.AA-27376
follow a particular scaling relation based on the Kolmogorov transportwithrespecttothebackgroundmagneticfieldB .
0
phenomenology given by k¯ (cid:96)1/3 = k¯2/3 (hereafter GS scaling)
(cid:107) ⊥
1. This relation relies on a critical balance between the perpen- This paper is organized as follows. Section 2 presents the
dicular cascade and Alfvén wave packets crossing timescales. generalnumericalframeworkadoptedinthispaper.Wedescribe
Inthecompressiblelimit,thesituationismorecomplex:Alfvén in§2.1theMHDsimulationsandtheforcingprocedureaswell
and slow-magnetosonic modes have been found to follow the as the set-up of our simulations. In §2.2 we describe the ki-
sameGSscaling,butfast-magnetosonicmodeshavebeenfound neticmoduledevelopedtointegratetheparticletrajectories.Sec-
tofollowaKraichnanspectrumandtohaveanisotropiccascade tion§3presentsthemainresultsobtainedinthispaper;in§3.3
(Cho&Lazarian2003). theCRmfpdependencewiththeAlfvénicMachnumberisde-
The impact of the different turbulence models over the CR scribed, while in §3.5 the dependence with particle rigidities is
transporthasbeendiscussedinalonglistofarticles.Apartfrom investigated.In§4wesummarizeandconcludeourwork.
the above-mentioned analytical calculations, another approach
based on numerical simulations uses synthetic turbulence
2. Kinetic-MHD simulations
modelsandthenreconstructsthedifferenttransportcoefficients
by averaging CR trajectories over a large number of particles In this work particle trajectories are calculated by means of di-
and magnetic cube realizations (Giacalone & Jokipii 1999; rect numerical simulation (DNS) of particle transport in MHD
Casse et al. 2002; Candia & Roulet 2004; Marcowith et al. turbulence. To this end we upgraded the RAMSES MHD code
2006; Tautz 2010; Sun 2011; Laitinen et al. 2012; Hussein (Teyssier2002;Fromangetal.2006)byincludingaturbulence
etal.2015).Inallcases,thesecalculationsimposeamodelfor forcingmodule(see§2.1)andakineticmodule(see§2.2).
the turbulent spectrum necessary to derive the particle mfps.
With the progress of computational power, some studies have
beenperformedtotestCRtrajectoriesinmorerealistic(orless 2.1. Magneto-hydrodynamicsimulations
model-dependent) situations where the turbulence spectrum
2.1.1. Forcingmodule
is calculated by directly solving some fundamental equations.
This is particularly the case in the context of CR acceleration The MHD simulations were performed in a 3D cartesian grid
at non-relativistic (and also relativistic) shocks where particle- withperiodicalboundaryconditions.TheHLLDRiemannsolver
in-cell (PIC) techniques calculate the electromagnetic field foridealMHDisusedinthiswork.Thesimulationboxhasasize
solutions from the Maxwell equations. Recent PIC simulations L=1thatcanberescaledafterwardtodefinethescalelengthsof
have started self-consistent investigations of particle injection theproblemunderconsideration.Theturbulenceisgeneratedby
andaccelerationintherelativisticregimeatsupernovaremnant forcingthefluidvelocitycomponentwithanexternalforcef in
shock fronts (Bai et al. 2014). At higher energies (and larger theEulerequation.TheforceisexpressedintermsofitsFourier
scales), calculations coupling MHD and PIC techniques have transform fˆ.Itisobtainedbytherealizationofindependentran-
permittedinsightsintothedynamicsoftheshockCRprecursor domOrnstein-Uhlenbeckprocessesforeachexcitedwavenum-
(Revilleetal.2008). ber.Adetaileddescriptionoftheimplementationoftheforcing
In this work, we also use a PIC-MHD approach to investigate isdescribedinSchmidtetal.(2009)butwereproduceitherefor
the propagation of CRs in magnetized turbulent media. The convenience.
turbulenceisproducedbydirectnumericalsimulationsusingthe The components of the force are expressed in terms of the
RAMSES MHD code. We follow a few pioneering calculations Fourier amplitudes of N modes (typically N = 32). One ad-
m m
performed in the context of space plasmas (Wisniewski et al. vantage of the Ornstein-Uhlenbeck process is that each mode
2012)orinterstellarmediumstudies(Beresnyaketal.2011;Xu component evolves statistically independently from the others.
&Yan2013).Wefirstsummarizetheresultsobtainedinthelat- Wehave
tertwoworksaswearemoreinvolvedinCRpropagationinthe
interstellarmedium.Wethenextendtheircalculationstoinclude fi =ΣmNm=1fˆi,mcos(2πkj,mxj), (1)
other turbulence regimes and discuss issues connected with the
turbulence forcing process. This work presents an investigation where fˆi,m areinitializedtozeroatthebeginningofeachMHD
of CR transport in MHD turbulence over a large parameter run.Eachmodefollowsastochasticdifferentialequation:
space of Alfvénic Mach numbers and particle rigidities. Our
 (cid:114) 
sciamseusl.atWioensn,ohtoewtehvaetr,inarethriesstwriocrtekdthtoesAubfv-AénlficvéMnicacthurbnuulmenbceer dfˆi,m =gχ×−(cid:15)¯fˆi,mdTt +β¯cTs 2wT(k)2Pχi,j,mdξj . (2)
M is defined as the ratio of the rms turbulent fluid velocity
a
V = (cid:104)δV2(cid:105)1/2 totheAlfvénspeedtakenwiththetotalmagnetic Thedifferentparametersenteringintheforcingarewtheampli-
(cid:113)
field V = B /(cid:112)4πρ with B = δB2+B2 and where ρ is tudeoftheforcing,T theauto-correlationtimescaleand(cid:15)¯andβ¯
a T T 0 areparameterscontrollingtherelativestrengthofthedetermin-
fluidmassdensity.Hencehereafter Ma < 1,butalsoδB/B0 ≤ 1 isticandstochasticforcingcomponents((cid:15)¯ =β¯ =1isadoptedin
so that we can always define a parallel and a perpendicular √
thiswork)anddξ = ξ dt.Therandomvariableξ issampled
j j j
over a Gaussian distribution with zero mean and variance one.
1 k¯andkaretheperturbedwavenumbersdefinedalongthelocalmag- Thetensor Pχi,j,m definesthegeometryoftheforcing(Federrath
neticfieldB andthetotalmagneticfieldB ,respectively.Wealsode- etal.2008,2011),
L T
fine the global or background magnetic B , i.e. the large scale back-
0 kk
groundmagneticfield.Thetotalmagneticfieldincludesthecontribu- Pχ =χδ +(1−2χ) i j . (3)
tionofthebackgroundmagneticfieldandtheperturbedmagneticfield i,j,m i,j k2
δB.Thelocalmagneticfieldisnotunique.Itsdirectionvarieswiththe
scale 1/k¯ under specification (see the discussion in Cho & Vishniac where χ ∈ [0,1] controls the relative importance of solenoidal
(cid:107)
(2000)).(cid:96)isthecoherencelengthoftheturbulence. andcompressivemodesintherandomforcingoperatorentering
Articlenumber,page2of15
R.Cohet andA.Marcowith:Cosmicraypropagationinsub-Alfvénicmagnetohydrodynamicturbulence
in the Ornstein-Uhlenbeck process. Hence, χ = 0 corresponds Name Level w χ (cid:104)η(cid:105) (cid:104)M (cid:105) (cid:104)M (cid:105)
s a
to a pure compressive forcing (hereafter CoF or curl-free) and 8J01c0.0 8 0.04 0.0 0.09 0.08 0.90
χ=1correspondstoapuresolenoidalforcing(hereafterSoFor 8J02c0.5 8 0.05 0.5 0.35 0.58 0.58
divergence-free). We also use a mixed forcing (hereafter MoF) 8J03c1.0 8 0.1 1.0 0.83 0.1 1.00
(cid:112)
withχ = 0.5.Thenormalizationfactorg = 3/ 1−2χ+3χ2. 9J01c0.0 9 0.04 0.0 0.25 0.35 0.34
χ
The fraction of compressible modes generated by a forcing 9J02c0.0 9 0.07 0.0 0.38 0.50 0.50
of strength χ in 3D is derived from the norm of each part of 9J03c0.0 9 0.09 0.0 0.50 0.61 0.66
the forcing tensor (right-hand side terms in Eq. 3). One gets 9J04c0.0 9 0.18 0.0 0.71 0.92 0.88
r = F /(F +F ) = (1−χ)2/(1−2χ+3χ2)(Federrath 9J05c0.5 9 0.02 0.5 0.18 0.37 0.37
c comp sol comp
etal.2010).Hencer =0,1/3,1forχ=1,0.5,0respectively. 9J06c0.5 9 0.045 0.5 0.34 0.58 0.58
c
MacLow&Klessen(2004)haveproposeddifferentsourcesof 9J07c0.5 9 0.08 0.5 0.52 0.68 0.67
interstellarturbulence:galacticspiraldensityshocks,large-scale 9J08c0.5 9 0.15 0.5 0.73 0.92 0.86
gravitational contraction, supernova explosions, and proto- 9J09c1.0 9 0.01 1.0 0.18 0.34 0.34
stellar jets and winds. All are likely to excite compressible 9J10c1.0 9 0.0225 1.0 0.32 0.50 0.50
modes. Incompressible forcing can, for instance be produced 9J11c1.0 9 0.045 1.0 0.49 0.69 0.66
via shearing flows. It is therefore not unrealistic to consider 9J12c1.0 9 0.1 1.0 0.72 1.0 0.88
that the turbulence forced at the largest scales has a mixture 10J01c0.5 10 0.045 0.5 0.35 0.53 0.53
of compressible and incompressible components. The effect of 10J02c1.0 10 0.01 1.0 0.18 0.37 0.37
compressibleandincompressibleforcingoverMHDturbulence
hasalsobeenaddressedinWisniewskietal.(2013). Table1:MHDsimulationsusedinthiswork
In all of our simulations the isothermal approximation for
η η
thegasequationofstateisused.Thepressureandthedensityin 0
0. 0.
the simulation set-up are P = ρ = 1 unless otherwise specified
√ 0.10 0.10
andthesoundspeedisc = γ (cid:39)1.005.Theenergyisinjected
s 0.30 0.31
intothelarge-scalemodesinthewavenumberintervalk∈]1,3]×
0.50 0.57
2π/N(whasaHeavisidefunctionshapeink),whereNisdefined
0.60 0.75
bythelevelofrefinementXasN =2X.Wenotethathereafterwe √
1/ 2∼0.707 1.00
makethedistinctionbetween L,whichisthesizeofsimulation
0.80 1.33
cube, and L < L, which is the scale of turbulence injection
inj
0.90 2.10
identifiedwiththecoherencelength(cid:96)oftheturbulence.
1.00 ∞
2.1.2. Simulationset-up Table2:Correspondencebetweentheparametersηandη .
0
We performed different types of MHD simulations; they are
summarizedintable1.Thenomenclaturefollowsfromthevalue
advantageofusingηisthatthisquantityisboundedbetween0
of the forcing parameter χ and the resolution level X. For in-
and1.Inourworkbothηand M areaveragedoverthreeMHD
stance,asimulationataresolution N = 5123 = 29 isatalevel a
snapshots and are written (cid:104)η(cid:105) ≤ (cid:104)M (cid:105). In table 1 the value of
X=9.Wethenspecifyχasc = 0,0.5,1and J thejobnumberto a
(cid:104)η(cid:105)forthedifferentjobsisgiveninthefifthcolumn.Comparing
differentiatejobswithidenticalχvaluesatagivenlevel.Weper-
thesevalueswiththosedispalyedintable2weseethatoursim-
formedsimulationsatlevels8,9and10.Alsodisplayedintable
ulationsareinthesub-totrans-Alfvénicregimewithrespectto
1aretheresultingmeansonicandAlfvénicMachnumbersob-
theAlfvénvelocitytakeninthebackgroundmagneticfield.
tainedaftertheturbulentspectrumhasreachedaquasi-stationary
state.ThevalueofthemeanAlfvénicandsonicMachnumbers
areaveragedoverthreesnapshotswhichwillbeusedasturbulent 2.2. Kineticmodule
fieldrealizationsforthepropagationofparticles.Eachsnapshot
isseparatedbyseveralcascadingtimescalesatlargescales(i.e. 2.2.1. Generaldescription
∼ 10L /V, where V is the perturbation velocity). The MHD
inj WesolvedtheLorentzequationforeachparticleofchargeqand
simulations have been performed at the CINES center on the
massm.Theparticlehasamomentump = γmvandanormal-
Jade and Occigen super-calculators. The typical duration of a
izedvelocityβ=v/candpropagatesinanelectromagneticfield
MHDrunat2563,5123,and10243 are5000,20000,and90000
δE(nomeanelectricfield),B =δB+B :
h.cpu,respectively. T 0
Beforepresentingourcalculationsindetailwewanttoclar-
dp
ifyanissue,barelyexplicited,concerningthedefinitionsofthe = qδE+qβ∧B ,
dt T
Alfvénic Mach number in the literature. In this work we use
Ma = V/Va with Va defined in the total magnetic field (see dr = v. (4)
§1) and V is the rms turbulent fluid velocity. Instead, several dt
other works use M = V/V > M with V defined in the
a0 a0 a a0
TheelectromagneticfluctuatingcomponentsδEandδBarepro-
background magnetic field. The turbulence is usually forced in
vided by the MHD code. In this work, only protons are con-
the velocity space so there is additional step to link M (or
a sidered, hence m = m and q = +e. Each particle is injected
M ) to the level of magnetic fluctuations. Here we consider p
a0 (cid:113) with a Lorentz factor γ . This defines the particle’s Larmor ra-
η = δB/ (δB2+B20)whereasinmanyotherworksη0 = δB/B0 diusrL0 =γ0mc2/eB0.W0 ealsodefinetheparticle’ssynchrotron
isused.Table2givesthecorrespondencebetweenηandη .The pulsation Ω = r /c. Eq. 4 is normalized with respect to this
0 0 L0
Articlenumber,page3of15
A&Aproofs:manuscriptno.AA-27376
initialvalueas(seeBeresnyaketal.(2011),hereafterBYL11): byBYL11andXu&Yan(2013)(hereafterXY13),bothdevel-
opedinthecontextofCRpropagationintheinterstellarmedium.
duˆ δE δB
= γˆ +uˆ∧( +e ), (5) In §3.3 we present our results describing CR mfp dependences
dtˆ B0 B0 z withrespecttoAlfvénicMachnumberswhilein§3.5wepresent
dxˆ r our results describing CR mfp dependences with respect to CR
= L0 ×uˆ . (6)
dtˆ L rigidities. In §3.4 we discuss how forcing geometries affect the
CRtransport.Finally§3.6presentsthevariationsoftheratioof
Weassumethatthelarge-scalemagneticfieldB isalignedalong
0 perpendiculartoparallelmfpswiththeparticlerigidityatdiffer-
e . In the previous equations we have the following notations:
z entturbulencelevels.
uˆ = γˆβ, γˆ = γ/γ , tˆ = t × (γˆr /c)−1. In practice, at least in
0 L0
theinterstellarmedium,theeffectoftheelectricfieldoverhigh-
energy(multiTeV)CRscanbeneglectedandγˆ remainsequalto 3.1. Summaryofpreviousworks
1.
Beresnyak et al. (2011) have investigated the propagation
When particles are submitted to radiative losses (or re-
of CRs in incompressible forced MHD turbulence. In the
acceleration) it is necessary to add another equation for γˆ. For
incompressiblelimit,theeffectsofthefast-magnetosonicmodes
instance,relativisticprotonscoolingundertheeffectofpionpro-
are absent, the transport is controlled by the shear-Alfvén and
ductionwithacrosssectionσ verify
pp
the pseudo-Alfvén modes (the incompressible limit of slow-
dlnγˆ r magnetosonicmodes).Bothtypesofmodeshavebeenfoundto
=−γˆ × L0 ,
dtˆ ct follow a GS scaling (Maron & Goldreich 2001). BYL11 solve
pp
Eqs. 5 and 6 neglecting the effect of electric field fluctuations
wheret (cid:39) (σ n c)−1 foramediumofhydrogendensityn . and use different statistically independent (static) magnetic
pp pp H H
Wenotethathereafter,astheparticleenergyisfixed,wedenote configurations. No back-reaction of CRs over the turbulence is
r andΩ tor andΩ. consideredinthiswork.
L0 0 L
Beresnyak et al. (2011) have conducted 3D simulations at a
resolution of 7683 in both balanced and imbalanced turbulence
2.2.2. Integrationschemes
cases. Imbalanced turbulence is expected, for instance close to
CRsourceswhereagradientofCRcandriveperturbationsina
We tested several integration schemes for Eqs. 5 and 6: a leap-
frog scheme (by default implemented in RAMSES to treat the preferreddirectionwithrespecttothemeanmagneticfieldorin
the solar wind. The turbulence is injected at large scales using
propagation of test particles in a gravitational field), a Runge-
a stochastic solenoidal forcing (see the definition in §2.1.1).
Kutta method of 5th order, and a Bulirsch-Stoer method (Press
The authors have carried out simulations in a regime of weak
et al. (1992)). We tested the different integration methods and
turbulence(or strongbackgroundfield) withanAlfvénic Mach
found differences at the level of a few percent over the parti-
number M =0.1andinaregimeofstrongturbulence(orweak
cle mfps. This small difference can be explained by the con- a
struction of the code: RAMSES adapts the time step in the PIC background field) with Ma = 1. We note that BYL11 calculate
the Alfvén velocity with respect to the background magnetic
moduleinawaythataparticlecannotcrossmorethanonegrid
field B , hence M has to be interpreted as M = δB/B
cell within one time step (see Teyssier (2002), §2.4). We used 0 a a0 0
in our work. The injection scale of the turbulence in their
theleap-frogsecond-orderintegrationmethodtosavecomputa-
simulation is (in box length units) L (cid:39) 0.2L. We note that in
tional resources. The essential effect now is the magnetic field inj
thesub-Alfvéniccase,BYL11useanelongatedboxwithamain
calculationattheparticleposition(seenextsection).Thetypical
durationofoneruntocompletethepropagationof106 particles axis parallel to B0. Then, in that case Linj corresponds to the
perpendicular turbulent outer-scale and the parallel outer-scale
depends on the particle’s rigidity and the Alfvénic Mach num-
ber. It varies from 106 to a few 107 time steps, in other words oftheturbulenceintheparalleldirectionis10Linj.
Beresnyak et al. (2011) have computed the variation of spatial
fromtenhourstoafewtensofhoursofcomputationontheJade
diffusion coefficients with respect to the direction of the back-
andOccigensupercomputersattheCINESfacility.
groundmagneticfieldB .Balancedandimbalancedturbulence
0
cases do not show strong differences and result in diffusion
2.2.3. Electromagneticfieldinterpolation coefficients varying by a factor of less than a few. In the case
of balanced trans-Alfvénic turbulence, BYL11 find a parallel
The magnetic field is interpolated at the position of the parti-
diffusion coefficient D independent of the particle’s rigidity
clefromitsvaluesderivedbytheMHDcodeatthegridpoints. (cid:107)
at low rigidity ρ = r /L ≤ 0.02 2 then scaling as ρ1/3 in the
For this, we have adopted a volume-averaged interpolation of L
range0.02 < ρ < 0.1andscalingasρ2 beyond(seetheirfigure
thefieldstrength.Thevolumecanincludethenexteightneigh-
6). Considering perpendicular diffusion, the authors claim to
bor points, this is the first-order cloud-in-cell (CIC) interpola-
tionimplementedbydefaultinRAMSESorthenext64neighbor have a coefficient independent of the particle’s rigidity in both
sub-Alfvénicandtrans-Alfvéniccases.Infact,consideringtheir
pointsusingathird-orderpiecewisecubicspline(PCS)interpo-
figures 7 and 8, their solutions show stronger dependence on
lation (Haugbølle et al. (2013)). The results in this work were
rigidity. Especially, D ∝ ρ1/3 fits the curve corresponding to
obtainedusingthePCSinterpolationmethod. ⊥
the balanced turbulence case in the trans-Alfvénic case better.
We find the same trend in the sub-Alfvénic case for ρ ≤ 0.05.
3. Direct numerical simulations of cosmic ray Noexpliciteffectsabouttheturbulencelevelhasbeenreported,
transport but a gain of a factor of ∼ 3 − 4 of difference for D⊥ can be
Direct numerical simulations of CR transport in magnetic tur- 2 Weadoptthefollowingnotations:ρ˜ isthereducedrigiditywithre-
bulence is a young but already extended field of research. In specttotheturbulenceouterscaleandρisthereducedrigidityincube
§3.1 we summarize the results of two related works proposed units;inBYL11andinourworkwehaveρ˜ ∼5ρ.
Articlenumber,page4of15
R.Cohet andA.Marcowith:Cosmicraypropagationinsub-Alfvénicmagnetohydrodynamicturbulence
deducedbetweensub-andtrans-Alfvéniccases.
Xu & Yan (2013) have investigated the propagation of
CRs in compressible MHD turbulence and so the simulations
retaintheeffectsofshear-Alfvénmodesandbothmagnetosonic
modes. Here again the magnetostatic and test-particle limits
havebeenretained.
Xu & Yan (2013) conducted 3D simulations at a resolution of
5123.Theturbulenceisinjectedatlargescalesusingastochastic
forcing of the velocity field. The forcing is purely solenoidal.
The authors carry simulations in a regime of weak turbulence
with Alfvénic Mach numbers ranging from M = 0.19 to
a
M = 0.73. XY13 have calculated the Alfvén velocity with
a
respect to the total magnetic field B hence their definition of
T
M is similar to our. The injection scale of the turbulence in
a
theirsimulationis(inboxlengthunits)L =0.4L.
inj
Xu & Yan (2013) computed spatial diffusion coefficients with Fig.1:Time-dependentevolutionofλ (t)(continuouslines)and
(cid:107)
respect to the direction of the background magnetic field B0. λ⊥(t)(dashedlines)foranormalizedparticleLarmorradiusρ=
The authors conducted a detailed study of the effect of the r /L=0.037atthreedifferentAlfvénicMachnumberswithχ=
L
turbulence level over the particle transport and have calculated 1(jobs9J09c1.0,9J11c1.0and9J12c1.0).Themfpcorresponds
CR parallel and perpendicular mfps dependences with respect to the plateau reached at large timescales. Time is in units of
to Ma. The results, however have been presented only at one particlegyroperiodτ=Ω−1.
particle rigidity ρ = 0.01 (so in box length units). Parallel
mfps have been found to be > L scaling roughly as M−2.
inj a
Perpendicularmfpshavebeenfoundtobeconsistentwitha M4 at ρ = r /L = 0.037). We note that at low M the ballistic
a L a
scaling expected theoretically by Yan & Lazarian (2008) in the regime is present up to t ∼ 103Ω−1. Once the diffusion coeffi-
casetheparallelmfpsverifyλ > L (seedetailsin§3.3.3). cient κ = lim κ (t) is found the corresponding mfp is ob-
(cid:107) inj (cid:107),⊥ t→+∞ (cid:107),⊥
tained using the relation λ = 3κ /v (the particle velocity
(cid:107),⊥ (cid:107),⊥
We present our results below. We provide an extensive in- is identified to the light velocity (v = c) in the cases consid-
vestigationofCRparallelandperpendicularmfpsdependences
ered here). For the results presented below a typical number of
withrespecttoMaandρthatextendthetwoaboveworkssignif- N =106particlesisused.
p
icantly.Wefindpower-lawvariationsforλ andλ .Thepower-
(cid:107) ⊥
lawindexisdenotedαhereafter.
3.3. CRmeanfreepaths:AlfvénicMachnumber
dependences
3.2. CRmeanfreepaths:Calculationprocedure
Wefirstpresentthedependenceofλ andλ withrespecttothe
(cid:107) ⊥
Cosmicraymfpsarereconstructedusinganaveragingprocedure
turbulencelevel M .Inthefollowingfigures,theselectedρare
a
overalargenumberofindividualparticletrajectories.Werecord
0.01,0.037,0.047,0.058,0.074,0.097,0.117,0.147, and 0.186
the position x,y,z in the 3D simulation cube for each particle: respectively. We note that all rigidities except for ρ = 0.01 can
z(t) determines the position of the particle with respect to the
beassociatedwithascaleintheinertialrangeoftheturbulence
backgroundmagneticfieldB , x(t),y(t)determinestheposition
0 ascanbeseeninfigure2.
oftheparticleinthedirectionsperpendiculartothebackground
Anotherimportantaspecttokeepinmindwhilecomparingthese
magneticfieldB .
0 results with the transport in the interstellar medium (ISM) or
Toproceedtomfpcalculations,wehaveselectedasetofN =3
r interplanetarymedium(IPM)isthatoursolutionsarerestricted
magnetic field realizations separated by at least two large-scale
to high-energy particles. To fix the ideas we consider typical
cascade times in order to ensure a statistical independence be-
magneticfieldsof5µGintheISMand50µGinthesolarwind.
tween two successive realizations. For each realization an en-
Hence,consideringturbulenceinjectionscalesexpressedinpar-
semble of Np particles are propagated until a convergence is secandastronomicalunits,wehaveρ˜ ∼0.21E B−1 L−1
obtained (see next). The parallel spatial diffusion coefficient is ISM PeV 5µG inj,pc
and ρ˜ ∼ 4.4 × 10−3E B−1 L−1 in the ISM and IPM,
calculatedusingthefollowingformula: IPM GeV 50µG inj,AU
respectively. Here E and E are the particle energies
PeV GeV
expressedinPeVandGeVunits.
κ (t)= (cid:104)(z(t)−z(0))2(cid:105) = 1 1 ΣNr ΣNp (zi,j(t)−zi,j(0))2 . (7)
(cid:107) 2t N N i=1 j=1 2t
r p
The perpendicular spatial diffusion coefficient is calculated 3.3.1. Parallelmeanfreepaths
usingthefollowingformula:
In figure 3, we show the dependence of the CR parallel mfps
with respect to the Alfvénic Mach number (the upper figure
κ⊥(t)= N1 N1 ΣiN=r1ΣNj=p1(xi,j(t)−xi,j(0))24+t(yi,j(t)−yi,j(0))2 . ddiissppllaayyss tthhee rreessuullttssoobbttaaiinneeddffoorrχχ ==10((SCooFF)),)t.hIetlcoawnebrefigseuerne
r p
that parallel mfps decrease with M at all particle rigidities.
(8) a
TheSoFsolutionsshowfastervariationswith M thantheCoF
a
These two coefficients are calculated until they both con- solutions. We find typical ratios λ (χ = 1)/λ (χ = 0) varying
(cid:107) (cid:107)
verge to a plateau (see an example in figure 1 for a particle between ∼ 100 at low M to ∼ 1 as M → 1. We find that
a a
Articlenumber,page5of15
A&Aproofs:manuscriptno.AA-27376
10000
1000
100 M-2
A
L
/λ| | MA-4
10
r =0.01
L
r =0.037
1 rL=0.074
L
r =0.097
L
r =0.147
L
0.1
1
M
a
100
r =0.01
L
r =0.037
L
r =0.074
L
Fig.2:Stationarykineticenergyturbulentspectracorresponding 10 rL=0.097
r =0.147
L
to jobs 8J01c0.5, 9J06c0.5, and 10J01c0.5 in table 1. Also dis-
played: the normalized wavenumbers kL corresponding to dif-
L
ferentnormalizedLarmorradiiρ=rL/L. /λ| | 1
λ⊥(χ = 1)/λ⊥(χ = 0) varies from 4-5 for Ma < 0.6 to ∼ 1 as 0.1 MA-2
M → 1. In figure 3 (the SoF case) it can also be seen that the
a M-4
parallel mfp decreases with an increasing rigidity, which is an A
unsualbehavior.Wecomebacktothisissuein§3.5.4.
0.01
1
We extracted a subset of solutions at two different Larmor
M
radius ρ = 0.01 and ρ = 0.097. The former allows us to make a
some comparisons with XY13, while the latter corresponds Fig.3:[Upperfigure]ParallelmfpversusAlfvénicMachnum-
to a particle Larmor radius in resonance with a wave number bers M atdifferentparticlerigiditiesinthecaseofasolenoidal
a
k0.097 deepintheinertialrangeoftheturbulence.Itcorresponds forcing (χ = 1). [Lower figure] Parallel mfp versus Alfvénic
to a scale reasonably far from the injection zone (typically Mach numbers M at different particle rigidities in the case of
a
k0.097L ∼ 10,k0.097Linj ∼ 2), as in our case Linj ∼ 0.2L. acompressibleforcing(χ = 0). M−2 and M−4 dependencesare
a a
Particles with rigidities in the inertial range are interesting for shownintheinsert.
investigating the effect of resonant pitch-angle scattering over
particleparallelmfps.Thedifficultyisthatresonantinteractions
occurringatkr ∼ 1areimportantinthecalculationofλ ,and
L (cid:107)
non-resonant interaction produced by the transit-time damping errorsobtainedusingthreeandsixsnapshotsaresimilarto(even
(TTD) mechanism due to fluctuations at scales larger than r if slightly reduced in the case with six snapshots), but are re-
L
also contribute (Lee & Voelk 1975). The error bars correspond ducedwithrespecttothecasewithonesnapshot.Inordertoop-
to standard deviations produced by the solutions from the timizetheaveragingprocedurewekeepthenumberofsnapshots
differentsnapshots. Nr =3inallourcalculations.
Inthecaseofincompressible(solenoidal)forcingχ = 1we
findthefollowingsolutions: 3.3.2. Perpendicularmeanfreepaths
λ(cid:107) ∝ Ma−7.69±0.57,ρ=0.01 The perpendicular mfps are controlled by two efects: a first
λ ∝ M−4.91±0.35,ρ=0.097 (9) contribution is associated with the different processes entering
(cid:107) a in the parallel transport (resonant and non-resonant pitch-angle
Inthecaseofcompressibleforcingχ = 0wefindthefollowing scattering) and a second contribution is associated with the
solutions: wanderingofmagneticfieldlines(Jokipii1971).Figure4shows
the CR perpendicular mfps as a function of the Alfvénic Mach
λ ∝ M−2.44±0.82,ρ=0.01
(cid:107) a number.
λ ∝ M−1.82±0.35,ρ=0.097 (10)
(cid:107) a
Inthecaseofsolenoidalforcingχ=1wefindthefollowing
Theerrorontheindexisobtainedconsideringtheextremeslopes
solutions:
that fit all data points. The mean index values and their errors
dependontheerrorsproducedbyeachdatapoint.Wehavecom-
pared our results using one, three and six snapshots for the job
λ ∝ M2.94±0.69,ρ=0.01
9J06c0.5. We found that mfps in each of the last two cases are ⊥ a
containedwithintheerrorofthesimulationatonesnapshot.The λ ∝ M2.52±0.26,ρ=0.097 (11)
⊥ a
Articlenumber,page6of15
R.Cohet andA.Marcowith:Cosmicraypropagationinsub-Alfvénicmagnetohydrodynamicturbulence
0.1
r =0.01 dominatestheCRperpendiculartransportitisnotunrealisticto
L
rL=0.037 obtainsimilarresultsforλ⊥inbothstudies.
rL=0.074 Atthisstageitisdifficulttoexplainthereasonforthediscrep-
0.01 rrL==00..019477 ancy found between the two works for λ(cid:107) at low and moderate
L M . What can be said is that even if the forcing is incompress-
a
ible in both cases, the forcing method is not similar in the
L
/ 0.001 two MHD codes. Also, the structure of the two MHD codes
⊥
λ and the methods used in the codes to derive fluid solutions are
probablydifferent.AnothercriticalaspectisthatXY13derived
0.0001 M4 their solutions at a normalized rigidity ρ = 0.01. Particles at
A
this rigidity have a Larmor radius that is in resonance with
M2 perturbations that are in the dissipation range of the turbulence
A
for simulations at a resolution of 5123 (as can be clearly seen
1e-05
1 in figure 2). Wave-particle resonance is a critical process that
controls particle scattering along the magnetic field lines. It
M
a can be seen especially in the SoF case (and to a lesser extent
0.1
rL=0.01 in the CoF case) as all rigidities except for ρ = 0.01 (and also
rrL==00..003774 ρ=0.037atlowMa)haveauniformdependenceonMa.Thisis
L themainreasonwhywithregardtotheproblemofpropagation
0.01 rL=0.097
r =0.147 ofCRinMHDturbulence-wesuggestthattheresults(hereand
L
inXY13)derivedatρ = 0.01havetobeinterpretedwithsome
/L⊥ 0.001 coathreerawtisleeasstpeccoinficeedr,niwngeλc(cid:107)o.nTcehnistriastealosourwahnyalhyesriseaoftner,Luanrmleossr
λ
radiithatcorrespondtoscalesintheinertialrangeoftheturbu-
lentspectrum(thatistosayallrigiditiesbutρ=0.01infigure2).
0.0001 M4
A
M2 Theoreticalframeworks:QLTpredictsλ(cid:107) ∝ Ma−2 becauseinthis
A approximation the scattering frequency is ∝ M2 (see e.g. Sun
a
1e-05 (2011))3.Casseetal.(2002)inthecaseofisotropicturbulence
1 have obtained λ ∝ η−2 ∝ M−2 in the limit δB (cid:28) B . They
(cid:107) a 0
Ma have found this scaling to be valid in the regime η → 1, i.e.
beyond the QLT validity domain (see their appendix A) 4 . A
Fig.4:[Upperfigure]PerpendicularmfpversusAlfvénicMach
deviation from the QLT result is also expected in the case of
numbers Ma at different particle rigidities in the case of a a slab-type turbulence with a steep inertial spectrum index 5 a
solenoidal forcing (χ = 1). [Lower figure] Perpendicular mfp
resultobtainedusingasecondorderquasi-lineartheory(Shalchi
versusAlfvénicMachnumbersM atdifferentparticlerigidities
a etal.2009).
inthecaseofacompressibleforcing(χ=0).
Forperpendicularmfps,QLTisobtainedinthelimitofvanish-
ing scattering (or infinite mfp) along the mean magnetic field
Inthecaseofcompressibleforcingχ = 0wefindthefollowing (see Jokipii & Parker (1968); Sun (2011)). In that regime, field
line wandering controls the perpendicular CR transport and
solutions:
leadstoλ ∝ M2.However,recenttheoreticalpredictionshave
⊥ a
λ⊥ ∝ Ma1.39±0.73,ρ=0.01 been proposed by Yan & Lazarian (2008) (see the discussion
λ ∝ M1.42±0.37,ρ=0.097 (12) in their §5) 6 in the framework of the compressible MHD
⊥ a turbulence model. Sub-Alfvénic turbulence is characterized by
tworegimes(Lazarian2006):betweenL and(cid:96) = L M2 the
inj tr inj a
turbulence is weak, the magnetic power spectrum develops in
3.3.3. Discussion
the perpendicular direction with respect to the mean magnetic
Comparison with previous works: We first compare our results field and scales as k−2; beyond (cid:96) up to scales such that
⊥ tr
intheSoFcase(χ=1)withthepublicationofXY13wherethe kL ∼ 40−50,whichfallinthedissipativerange,GSscaling
inj
mfp was calculated at ρ = 0.01 only. Comparing our figure 3 appliesandtheturbulenceisstrong.TheCRperpendicularmfp
(upper figure, red points) with figure 5 in XY13 we find some is mostly controlled by the wandering of magnetic field lines
discrepanciesespeciallyintheregime Ma < 0.7wherewefind but modulated by the effect of scattering: if λ(cid:107)/Linj > 1, then
λ /L > 100, whereas XY13 has typical mfps λ /L < 10. At perturbations uncorrelated at scales larger than (cid:96) produce a
(cid:107) (cid:107) tr
higher Alfvénic Mach numbers the two results are compatible. perpendicularmfpλ⊥/Linj = 1/3Ma4.Ifconverselyλ(cid:107)/Linj < 1,
Thiseffectisduetothemuchlargerindexα∼−7.7obtainedin then diffusion along the field lines lead to λ = λ M4. Finally
⊥ (cid:107) a
our simulations (see Eq. 9). Our CoF results at ρ = 0.01 are in we note that Casse et al. (2002) in the case of isotropic Kol-
goodagreementswithXY13solutions. mogorovturbulencenumericallyobtainedλ ∝ η4.6λ ,i.e.from
⊥ (cid:107)
Ifwenowcomparetheperpendicularmfpresults(theredpoints
inourupperfigure4andfigure6inXY13)wefindreasonably
3 Jokipii&Parker(1968)insteadpredictedλ ∝η−2andλ ∝η2.
goodagreement:ourindexisα=2.94±0.69andiscompatible (cid:107) 0 ⊥ 0
4 Casseetal.(2002)definedηasaquantitywhichcorrespondstoη2in
withtheXY13indexα = 4.21±0.75.Thenormalizationfactor
ourtext.
of the mfps is in agreement since we normalize λ⊥ to the box 5 Thisisthe1Dwavenumberindex,denoted sin§3.5.4.Deviations
length L, whereas XY13 normalized λ⊥ to the injection scale fromQLTareobtainedfors>2.
of the turbulence Linj. If the wandering of magnetic field lines 6 Yan&Lazarian(2008)usedMa0insteadofMa.
Articlenumber,page7of15
A&Aproofs:manuscriptno.AA-27376
theabovediscussionλ ∝η2.6. Second,iftheSoF-to-CoFmfpratioisreducedbyafactorof∼
⊥
Lazarian & Yan (2014) explored the effect of chaotic magnetic 10, then the gyro-resonance process by Alfvén and fast waves
field behavior over CR transport. The authors confirmed that contributes substantially. At higher Alfvénic Mach numbers
the CR perpendicular transport is diffusive in the regime the differences between the two forcing solutions vanish (the
λ /L > 1 for curvilinear distances s along the magnetic field ratioλ (χ = 1)/λ (χ = 0) → 1whatevertherigidityregime).It
(cid:107) inj (cid:107) (cid:107)
lines larger than a few tens of L . However, they described is difficult at this stage to draw conclusion from this effect. It
inj
the existence of a super-diffusive perpendicular transport could be associated with a stronger mode conversion between
regime produced by the strong turbulence regime at intermedi- Alfvén and fast-magnetosonic waves at scales smaller than the
aryscaless/L ∼1−10(seeEq.10inLazarian&Yan(2014)). injectionscales.
inj
However, the above analysis is still limited and a quantitative
Parallelmeanfreepath:Fromtheupperpartoffigure3itcanbe analysis would require to isolate the effects of the different
seen that our SoF results are incompatible with the predictions modes that compose the turbulent spectrum in the different
of QLT as α < −2. Even if we remove the two highest M forcing limits. One possibility would be to use the mode de-
a
points where QLT is questionable, we still find α < −2. We compositionprocedureproposedbyCho&Lazarian(2003)and
will consider this aspect again in §3.5.4 while discussing the improved by Kowal & Lazarian (2010), which should be valid
magnetic power spectra in the SoF limit. In the CoF case, we inthelimit Ma < 1(seealsothecasesdiscussedinWisniewski
obtain parallel mfps with indices close to -2 and this even etal.(2012)).
including the highest M points. This means that our results
a
are consistent with QLT predictions in the low M regime and
a
extenditinthehigh M regime.OurCoFresultsareconsistent Perpendicular mean free paths: Considering now figure 4, we
a
withthetheoreticalpredictionsofCasseetal.(2002). can first compare our results with the solutions obtained by
Comparing the upper and the lower parts of figure 3, it is clear Yan & Lazarian (2008). In the SoF case, λ(cid:107)/Linj varies in the
that SoF (χ = 1 ) simulations produce a turbulent spectrum range 1−103 (see the upper part of figure 3). At ρ = 0.01 we
that leads to a less efficient CR scattering (see also §3.5.4). have λ(cid:107)/Linj > 100 and λ⊥ ∝ Ma2.94±0.62. So the result is not
TheSoFfavorstheproductionofincompressibleAlfvénmodes compatible with a M4 scaling, but is still close to it. At such
a
with respect to compressible (in particular fast-magnetosonic) low rigidities and because λ /L (cid:29) 1, particles are expected
(cid:107) inj
modesevenifbothtypesofmodesaredrivenatlargescales.At to move along a field line over several coherence lengths with
smallerscales,theproductionoffastmodesfromAlfvénmodes a curvilinear abscissa s ∼ vt and their diffusion perpendicular
has been found to be weak at low levels of turbulence (Cho & mean magnetic field is determined by the divergence of mag-
Lazarian2003).OnecanexpectthattheSoFsolutionsatlowM neticfieldlines.Again,becausewave-particleinteractionsoccur
a
provide a good hint of the effect of Alfvén waves over the CR in the turbulence dissipation range at these rigidities, it is not
transport. Alfvén waves at the lowest perturbation order do not possible to fully trust that the derived parallel mfp is the one
contribute to TTD but only to Landau-synchrotron resonance to be expected if the turbulence were fully developed at these
(alsoknownasgyro-resonance),soparallelmfpsintheSoFcase scales. Including all the missing modes would likely produce a
likely result from this process. Recent transport models (see smaller parallel mfp. In turn, we cannot discard the possibility
e.g. Chandran (2000)) concluded that Alfvénic turbulence does that particle scattering could modify the result obtained for λ
⊥
not produce an efficient CR scattering because of the strong substantiallyevenifweexpecttheeffecttobesmallatlow M .
a
anisotropy the perturbations develop towards small scales. At We do not find any clear trend of super-diffusive CR transport
thesescales,theresonantconditionimpliesthatwithinaLarmor which could produce a larger final perpendicular mfps with
radiusaparticleinteractswithseveraluncorrelatedperturbations respect to the diffusive regime. Super-diffusion is produced by
resulting in an inefficient pitch-angle scattering. As reported the strong-Alfvénic turbulence (Lazarian & Yan 2014) in an
above,wefindtypicalratiosλ (χ = 1)/λ (χ = 0) ∼ 10−100at intervalwhichwidensas M increaseswhereasperturbedscales
(cid:107) (cid:107) a
low M depending on the particle rigidity. This partly supports associated with weak-Alfvénic turbulence which produces
a
the conclusion of the above analysis, but in the mean time the magnetic field line diffusion are restricted to scales ∼ L
inj
ratioisnotassmallaspredictedbythismodel.Yan&Lazarian only. We think that we need a higher grid resolution at low
(2002) have improved Chandran’s model by considering an M in order to include scales where GS spectrum develops
a
Alfvénic magnetic tensor derived from numerical calculations in order to capture a possible effect of super-diffusion. At
of Cho & Lazarian (2003) and found a parallel mfp larger by ρ = 0.097, we find λ ∝ M2.52±0.26, so we have a slower
⊥ a
∼4ordersofmagnitudeabove∼ 1GeV,whichisclosertoour Alfvénic Mach dependence at this rigidity. Again, the result is
results. It is also possible that a fraction of the energy of the not compatible with a M4 dependence. At high rigidities there
a
drivenvelocityfieldistransferredintofastmodesand/orthatthe is a competing effect between magnetic field line wandering
gyro-resonanceeffectinAfvénicturbulenceisstrongerthanex- and particle magnetic scattering. We have λ /L ∼ 1, and not
(cid:107) inj
pected.Forρ < 0.03,theparticlesarescatteredmainlybecause > 1, so wave-particle interactions have stronger effect over the
ofthemirroringeffectproducedbythelargestperturbedwave- perpendicular transport. This may explains why our solutions
lengths.ThedifferencesbetweenSoFandCoFsolutionscanbe arenotcompletelycompatiblewithλ ∝ M4.Itisagaindifficult
⊥ a
used to qualitatively estimate the impact of TTD over the CR to isolate an effect of super-diffusion in that rigidity regime.
transport:inaCoFgeometry,thefractionoffast-magnetosonic At high M , we find a perpendicular mean free path about two
a
waves is larger and the TTD contributes more to the scattering times smaller than the prediction λ /L ∼ 1/3M4 in relative
⊥ inj a
of particles. At higher rigidities, i.e. those interacting with agreement with the limit found by Yan & Lazarian (2008) and
turbulentperturbationsintheinertialrange,wecanadvancetwo in agreement with the results presented in XY13. In the CoF
arguments.First,consideringonlySoFsolutionsweseethatthe case, our solution is λ ∝ M1.42±0.37 (lower part of figure 4).
⊥ a
parallel mfp is reduced by a factor of ∼ 10 (with respect to the As λ is reduced with respect to the SoF case we also tested
(cid:107)
low rigidity regime); this shows that gyro-resonance by Alfvén the relation λ = λ M4. We find that our result is also not
⊥ (cid:107) a
wavescancontributetothepitch-anglescatteringsubstantially. completely compatible with the latter relation (although not
Articlenumber,page8of15
R.Cohet andA.Marcowith:Cosmicraypropagationinsub-Alfvénicmagnetohydrodynamicturbulence
1000
so far because the extreme values of α are 1.79 and 1.85),
respectively,butitisstronglyincompatiblewithλ ∝ M4.Here
⊥ a
again,theparallelmfpscoverarange(∼ 1−50)L .Thismay
inj
explain why we do not find a clear trend. The CoF cases are 100
characterizedbyastrongereffectofparticlescatteringproduced
byfast-magnetosonicwaves.
L
In both forcing geometries our solutions bracket the QLT / 10
solution λ ∝ M2 (from above for the SoF solutions and from λ| |
⊥ a
below for the CoF solutions). Our results are not completely
compatible with such a scaling possibly because the last point 1
at Ma = 0.89 is a bit under (above) the extrapolation of Ma2 in χ=0
the CoF case (SoF case). At high M indeed the use of QLT is χ=0.5
a
questionable. χ=1
0.1
1
Summary: Our solutions are globally consistent with the re- M
a
sults obtained by XY13 to the notable exception of λ(cid:107) at Ma < 0.1
0.7. This effect is possibly related to differences in forcing χ=0
χ=0.5
procedures and methods used to solve the MHD equations in
χ=1
the codes. If we accept that SoF (CoF) drives at large scales
Alfvén(fast-magnetosonic)wavespreferentially,ourresultscan
give some qualitative hints about the respective effect of gyro-
resonance and TTD processes. At low rigidity, i.e. for particles /L 0.01
inresonancewithmodesindissipationrangeoftheturbulence, ⊥
λ
the TTD likely controls the parallel CR mfp. However at such
rigidities the solutions have to be taken with care as resonant
modes are missing at small scales owing to the finite resolu-
tionofthesimulations.Withrespecttothelowrigidityfindings,
at higher rigidities the mfp is considerably reduced in the SoF
case, which can be due to the effect of gyro-resonance by the 0.001
1
Alfvénicturbulence.Also,theratiooftheparallelSoFandCoF
mfpsdecreases,whichcanbeduetotheeffectofgyro-resonance Ma
of Alfvén and fast-magnetosonic waves over the CR transport.
Wefindthatλ ∝ M−2 predictedbytheQLTiscompatiblewith Fig.5:[Upperfigure]:Parallelmeanfreepathdependencewith
(cid:107) a M at ρ = 0.097 for χ = 0,0.5,1. [Lower figure]: Perpen-
CoFsolutions,butnotwithwithSoFsolutions.Ourresultscon- a
dicular mean free path dependence with M at ρ = 0.097 for
cerningλ arelessaffirmativethantheonesadvancedbyXY13 a
⊥ χ=0,0.5,1.
and do not show a clear trend at the moment. The results have
not been found to be fully compatible with the theoretical pre-
dictions proposed in Yan & Lazarian (2008). One explanation
is that for each specific χ value the parallel mfp is not entirely
obtainedforχ=0.Thisfitswiththeargumentsadvancedabove
in the regime λ /L < 1 or λ /L > 1. Finally, our solutions
(cid:107) inj (cid:107) inj thatAlfvénwavesarepreferentiallyproducedwithrespecttofast
arebetweenthetwoscalings:λ ∝ M4 andλ ∝ M2 (thelatter
⊥ a ⊥ a magnetosonicwavesintheSoFandMoFgeometriesandthata
expectedfromQLT).Wedidnotfindacleartrendinourresults
tinyfractionoftheAlfvénenergyistransferredtofastwavesat
that could be associated with an effect of super-diffusion espe-
lowturbulentlevels(Cho&Lazarian2003).AthigherAlfvénic
ciallyatlowrigiditiesandlowAlfvénicMachnumberslikelybe-
numberstheeffectoftheforcinggeometryislessstringentpos-
causeGSturbulenceisrestrictedtoreducedscalelengths.More
sibly because of a more efficient mode conversion. There, we
simulationswithlargerdynamicsareneededinordertochoose
findthatparallelandperpendicularmfpsareindependentofχ.
definitivelyfromamongtheeffectsofparticlescatteringandany
effectofmagneticfieldlinetransportandalsotodecideamong
thedifferenttheoreticalmodels.Inparticular,resultsatM ∼0.1
a 3.5. CRmeanfreepaths:Rigiditydependences
wouldgreatlyhelptoprobetheQLTpredictionsfurther.Unfor-
tunately,suchlowAfvénicMachnumbersrequireverylongin-
3.5.1. Resolutiontests
tegrationtimescalestoreachadiffusiveplateau.Wecomeback
tothisaspectin§3.5.4.
Before discussing any specific rigidity dependence we first test
ourcalculationprocedureforρ > 0.037atdifferentgridresolu-
3.4. CRmeanfreepaths:forcingeffects tion. Figure 6 presents the parallel and perpendicular mfp with
respect to ρ at X= 8, 9, and 10. All the simulations were per-
Figure 5 presents the dependence of parallel and perpendicular formed at χ = 0.5 and jobs 8J02c0.5, 9J06c0.5, and 10J01c0.5
mfpsonMaatρ=0.097forχ=0,0.5,and1. wereused.
Itcanbeseenthatinbothcases(parallelandperpendicular Itcanbeseenthattheresultsarecompatibleeachotherwhatever
mfps),especiallyforM ≤0.5,theMoFsolutionsatχ=0.5are thevalueofX.Wealsoseeatrendwithadecreasingparallelmfp
a
closertotheSoFsolutions.Thisissupportedbythevalueofr astheresolutionincreases,especiallyatrigiditiesρ<0.1.Anef-
c
andthefractionofcompressiblemodes(see§2.1.1):intheMoF fectthatisexpectedsinceasXincreasestheinertialturbulence
caser = 1/3isclosertor = 0obtainedforχ = 1thanr = 1 zoneisshiftedtowardssmallerscales.
c c c
Articlenumber,page9of15
A&Aproofs:manuscriptno.AA-27376
100 10000
L=10 M =0.34
a
L=9 M =0.50
a
L=8 M =0.66
1000 a
1 1/2 Ma=0.88
1/3
0
100
L L -1
/ 10 /
λ| | λ| |
10
1
1 0.1
0.001 0.01 0.1 1 0.001 0.01 0.1 1
r /L r /L
L L
0.01 100
10
L L
λ /⊥ 0.001 /λ| |
1 1
1/2
M =0.34
a
M =0.50
L=10 1/3 Ma=0.66
L=9 a
L=8 0 Ma=0.88
0.0001 0.1
0.001 0.01 0.1 1 0.001 0.01 0.1 1
r /L r /L
L L
Fig.6:[Upperfigure]:Parallelmfpsatthreedifferentresolution Fig. 7: [Upper figure] Parallel mfps versus particle normalized
forjobs8J02c0.5,9J06c0.5,and10J01c0.5.[Lowerfigure]:The Larmorradiiρatdifferent M intheSoFcase(χ = 1).[Insets]
a
samebutforperpendicularmfps. Slopesobtainedfordifferentindices.[Lowerfigure]Parallelmfp
versusparticlenormalizedLarmorradiusρatdifferentM inthe
a
CoFcase(χ=0).
3.5.2. Parallelmeanfreepaths
3.5.3. Perpendicularmeanfreepaths
Figure7presentsCRparallelmfpsdependenceswithrespectto
theparticlenormalizedLarmorradiusρ=rL/L. Figure8presentstheCRperpendicularmfpwithrespecttothe
particlenormalizedLarmorradiusρ=r /L.
L
Intable3wepresentthefitindicesandtheerrorscorrespond-
ing to linear fits of the parallel mfp for 0.01 < ρ < 0.1 and Intable4wepresentthefitindicesandtheerrorscorrespond-
0.03<ρ<0.1inbothSoFandCoFcasesrespectively. ingtolinearfitsoftheperpendicularmfpfor0.01<ρ<0.1and
0.03<ρ<0.1intheSoFandCoFcasesrespectively.
Table 3: Parallel mean free path indices anderrors for SoF and
Table4:Perpendicularmeanfreepathindicesanderrorsforboth
CoFgeometries.
forcinggeometries.
Jobno. M α(0.03−0.1)±∆α α(0.01−0.1)±∆α
a Jobno. M α(0.03−0.1)±∆α α(0.01−0.1)±∆α
8J02c1.0 1.00 0.53±0.35 0.34±0.16 a
8J02c1.0 1.00 0.56±0.36 0.59±0.19
9J09c1.0 0.34 −1.19±0.36 −0.98±0.15
9J09c1.0 0.35 0.73±0.35 0.78±0.24
9J10c1.0 0.50 −0.21±0.35 −0.36±0.20
9J10c1.0 0.50 0.76±0.35 0.76±0.18
9J11c1.0 0.67 0.27±0.47 0.02±0.22
9J11c1.0 0.67 0.67±0.36 0.68±0.21
9J12c1.0 0.89 0.27±0.39 0.13±0.23
9J12c1.0 0.89 0.61±0.36 0.59±0.25
10J02c1.0 0.37 −0.87±0.35 −0.64±0.13
10J02c1.0 0.37 0.79±0.34 0.78±0.23
9J01c0.0 0.35 0.34±0.39 0.16±0.20
9J01c0.0 0.34 0.70±0.34 0.86±0.19
9J02c0.0 0.50 0.44±0.37 0.29±0.21
9J02c0.0 0.50 0.73±0.34 0.86±0.21
9J03c0.0 0.66 0.54±0.35 0.37±0.19
9J03c0.0 0.66 0.85±0.35 0.94±0.18
9J04c0.0 0.89 0.62±0.37 0.44±0.18
9J04c0.0 0.89 0.82±0.35 0.89±0.22
Articlenumber,page10of15