Cosmic ray propagation in subAlfvénic magnetohydrodynamic turbulence
Laboratoire Univers et particules de Montpellier, Université de
Montpellier, CNRS/IN2P3,
CC 72, Place Eugène Bataillon,
34095
Montpellier Cedex 5,
France
email:
Alexandre.Marcowith@umontpellier.fr
Received: 16 September 2015
Accepted: 13 January 2016
Context. The propagation of cosmic rays or energetic charged particles in magnetized turbulence is a complex problem which involves nonlinear waveparticle interactions and chaotic magnetic field lines transport. This problem has been addressed until recently using either analytical calculations or simulations using prescribed turbulence models. With the advent of super computers it is now possible to investigate energetic charged particle propagation using direct computation of electromagnetic fields. This is in particular the case for highenergy particles propagation in magnetohydrodynamic turbulence.
Aims. This work has the main objective to provide a detailed investigation of cosmic ray propagation in magnetohydrodynamic turbulent fields generated by forcing the fluid velocity field at large scales. It provides a derivation of the particle mean free path dependences in terms of the turbulence level described by the Alfvénic Mach number and in terms of the particle rigidity.
Methods. We use an upgrade version of the magnetohydrodynamic code RAMSES which includes a forcing module and a kinetic module and solve the Lorentz equation for each particle. The simulations are performed using a 3 dimension periodical box in the testparticle and magnetostatic limits. The forcing module is implemented using an OrnsteinUhlenbeck process. An ensemble average over a large number of particle trajectories is applied to reconstruct the particle mean free paths.
Results. We derive the cosmic ray mean free paths in terms of the Alfvénic Mach numbers and particle reduced rigidities in different turbulence forcing geometries. The reduced particle rigidity is ρ = r_{L}/L where r_{L} is the particle Larmor radius and L is the simulation box length related to the turbulence coherence or injection scale L_{inj} by L ~ 5 L_{inj}. We have investigated with a special attention compressible and solenoidal forcing geometries.
Conclusions. We find that compressible forcing solutions are compatible with the quasilinear theory or more advanced nonlinear theories which predict a rigidity dependence as ρ^{1/2} or ρ^{1/3}. Solenoidal forcing solutions at least at low or moderate Alfvénic numbers are not compatible with the above theoretical expectations and require more refined arguments to be interpreted. It appears especially for Alfvénic Mach numbers close to one that the wandering of field lines controls perpendicular mean free path solutions whatever the forcing geometry.
Key words: magnetohydrodynamics (MHD) / turbulence / methods: numerical / cosmic rays
© ESO, 2016
1. Introduction
Cosmic rays (CRs) propagate in our Galaxy through the interaction with electromagnetic perturbations usually described in the magnetohydrodynamic (MHD) approximation (Schlickeiser 2002), i.e. longwavelength perturbations with scales comparable to the particle’s Larmor radius r_{L} = vγ/ Ω_{c}, where Ω_{c} = ZeB/mc is the cyclotron pulsation of a particle of velocity v, mass m, and charge Z in a magnetic field of strength B. The modeling of particle transport in MHD turbulence is complex. It requires approximations to obtain the solutions that describe the particle trajectories. Transport studies were originally based on the quasilinear theory (QLT; Jokipii 1966; Schlickeiser 2002) where the unperturbed (gyromotion) particle trajectory is retained to derive the electromagnetic correlation tensor, itself used to derive the random Lorentz force exerted on the particles. The QLT is applicable over restricted timescale ranges, i.e. intermediate between pitchangle particle scattering and particle distribution isotropization times. The main drawback of the QLT is that it leads to infinite cosine pitchangle (the angle between the particle velocity and the background magnetic field) diffusion coefficients at 90 deg and hence produces a pathological behavior in the calculation of the particle parallel mean free path (mfp). Several analytical attempts have been proposed to deal with this issue (see Shalchi 2009 and references therein). In particular, the 90 deg scattering problem has been treated by the mean of the broadening of the resonance (Völk 1973). However, the way the resonance is broadened and the nonlinear solutions which have to be retained in diffusion coefficient calculations, are model dependent. Needless to say that all these studies are performed in the testparticle approximation (see Ptuskin et al. 2006). It is known that if the density in the relativistic particles is high enough, the particles can generate their own waves and induce a modification of the turbulence (Farmer & Goldreich 2004; Yan & Lazarian 2008, 2011).
Another difficulty in CR transport studies lies in the description of the random Lorentz forces. In the incompressible MHD limit the turbulence can be described by an anisotropic model (Goldreich & Sridhar 1995; Maron & Goldreich 2001). In this model, owing to the Alfvén wave packet dynamics, the turbulent perturbations are elongated along the magnetic field and follow a particular scaling relation based on the Kolmogorov phenomenology given by (hereafter GS scaling)^{1}. This relation relies on a critical balance between the perpendicular cascade and Alfvén wave packets crossing timescales. In the compressible limit, the situation is more complex: Alfvén and slowmagnetosonic modes have been found to follow the same GS scaling, but fastmagnetosonic modes have been found to follow a Kraichnan spectrum and to have an isotropic cascade (Cho & Lazarian 2003).
The impact of the different turbulence models over the CR transport has been discussed in a long list of articles. Apart from the abovementioned analytical calculations, another approach based on numerical simulations uses synthetic turbulence models and then reconstructs the different transport coefficients by averaging CR trajectories over a large number of particles and magnetic cube realizations (Giacalone & Jokipii 1999; Casse et al. 2002; Candia & Roulet 2004; Marcowith et al. 2006; Tautz 2010; Sun 2011; Laitinen et al. 2012; Hussein et al. 2015). In all cases, these calculations impose a model for the turbulent spectrum necessary to derive the particle mfps. With the progress of computational power, some studies have been performed to test CR trajectories in more realistic (or less modeldependent) situations where the turbulence spectrum is calculated by directly solving some fundamental equations. This is particularly the case in the context of CR acceleration at nonrelativistic (and also relativistic) shocks where particleincell (PIC) techniques calculate the electromagnetic field solutions from the Maxwell equations. Recent PIC simulations have started selfconsistent investigations of particle injection and acceleration in the relativistic regime at supernova remnant shock fronts (Bai et al. 2015). At higher energies (and larger scales), calculations coupling MHD and PIC techniques have permitted insights into the dynamics of the shock CR precursor (Reville et al. 2008).
In this work, we also use a PICMHD approach to investigate the propagation of CRs in magnetized turbulent media. The turbulence is produced by direct numerical simulations using the RAMSES MHD code. We follow a few pioneering calculations performed in the context of space plasmas (Wisniewski et al. 2012) or interstellar medium studies (Beresnyak et al. 2011; Xu & Yan 2013). We first summarize the results obtained in the latter two works as we are more involved in CR propagation in the interstellar medium. We then extend their calculations to include other turbulence regimes and discuss issues connected with the turbulence forcing process. This work presents an investigation of CR transport in MHD turbulence over a large parameter space of Alfvénic Mach numbers and particle rigidities. Our simulations, however, are restricted to subAlfvénic turbulence cases. We note that in this work the Afvénic Mach number M_{a} is defined as the ratio of the rms turbulent fluid velocity V = ⟨ δV^{2} ⟩ ^{1/2} to the Alfvén speed taken with the total magnetic field with and where ρ is fluid mass density. Hence hereafter M_{a}< 1, but also δB/B_{0} ≤ 1 so that we can always define a parallel and a perpendicular transport with respect to the background magnetic field B_{0}.
This paper is organized as follows. Section 2 presents the general numerical framework adopted in this paper. We describe in Sect. 2.1 the MHD simulations and the forcing procedure as well as the setup of our simulations. In Sect. 2.2 we describe the kinetic module developed to integrate the particle trajectories. Section 3 presents the main results obtained in this paper; in Sect. 3.3 the CR mfp dependence with the Alfvénic Mach number is described, while in Sect. 3.5 the dependence with particle rigidities is investigated. In Sect. 4 we summarize and conclude our work.
2. KineticMHD simulations
In this work particle trajectories are calculated by means of direct numerical simulation (DNS) of particle transport in MHD turbulence. To this end we upgraded the RAMSES MHD code (Teyssier 2002; Fromang et al. 2006) by including a turbulence forcing module (see Sect. 2.1) and a kinetic module (see Sect. 2.2).
2.1. Magnetohydrodynamic simulations
2.1.1. Forcing module
The MHD simulations were performed in a 3D cartesian grid with periodical boundary conditions. The HLLD Riemann solver for ideal MHD is used in this work. The simulation box has a size L = 1 that can be rescaled afterward to define the scale lengths of the problem under consideration. The turbulence is generated by forcing the fluid velocity component with an external force f in the Euler equation. The force is expressed in terms of its Fourier transform . It is obtained by the realization of independent random OrnsteinUhlenbeck processes for each excited wave number. A detailed description of the implementation of the forcing is described in Schmidt et al. (2009) but we reproduce it here for convenience.
The components of the force are expressed in terms of the Fourier amplitudes of N_{m} modes (typically N_{m} = 32). One advantage of the OrnsteinUhlenbeck process is that each mode component evolves statistically independently from the others. We have (1)where are initialized to zero at the beginning of each MHD run. Each mode follows a stochastic differential equation: (2)The different parameters entering in the forcing are w the amplitude of the forcing, T the autocorrelation timescale and and are parameters controlling the relative strength of the deterministic and stochastic forcing components ( is adopted in this work) and d. The random variable ξ_{j} is sampled over a Gaussian distribution with zero mean and variance one. The tensor defines the geometry of the forcing (Federrath et al. 2008, 2011), (3)where χ ∈ [ 0,1 ] controls the relative importance of solenoidal and compressive modes in the random forcing operator entering in the OrnsteinUhlenbeck process. Hence, χ = 0 corresponds to a pure compressive forcing (hereafter CoF or curlfree) and χ = 1 corresponds to a pure solenoidal forcing (hereafter SoF or divergencefree). We also use a mixed forcing (hereafter MoF) with χ = 0.5. The normalization factor . The fraction of compressible modes generated by a forcing of strength χ in 3D is derived from the norm of each part of the forcing tensor (righthand side terms in Eq. (3)). One gets r_{c} = F_{comp}/ (F_{sol} + F_{comp}) = (1 − χ)^{2}/ (1−2χ + 3χ^{2}) (Federrath et al. 2010). Hence r_{c} = 0,1/3,1 for χ = 1,0.5,0 respectively.
Mac Low & Klessen (2004) have proposed different sources of interstellar turbulence: galactic spiral density shocks, largescale gravitational contraction, supernova explosions, and protostellar jets and winds. All are likely to excite compressible modes. Incompressible forcing can, for instance be produced via shearing flows. It is therefore not unrealistic to consider that the turbulence forced at the largest scales has a mixture of compressible and incompressible components. The effect of compressible and incompressible forcing over MHD turbulence has also been addressed in Wisniewski et al. (2013).
In all of our simulations the isothermal approximation for the gas equation of state is used. The pressure and the density in the simulation setup are P = ρ = 1 unless otherwise specified and the sound speed is . The energy is injected into the largescale modes in the wave number interval k ∈ ] 1,3 ] × 2π/N (w has a Heaviside function shape in k), where N is defined by the level of refinement X as N = 2^{X}. We note that hereafter we make the distinction between L, which is the size of simulation cube, and L_{inj}<L, which is the scale of turbulence injection identified with the coherence length ℓ of the turbulence.
2.1.2. Simulation setup
We performed different types of MHD simulations; they are summarized in Table 1. The nomenclature follows from the value of the forcing parameter χ and the resolution level X. For instance, a simulation at a resolution N = 512^{3} = 2^{9} is at a level X = 9. We then specify χ as c = 0,0.5,1 and J the job number to differentiate jobs with identical χ values at a given level. We performed simulations at levels 8, 9 and 10. Also displayed in Table 1 are the resulting mean sonic and Alfvénic Mach numbers obtained after the turbulent spectrum has reached a quasistationary state. The value of the mean Alfvénic and sonic Mach numbers are averaged over three snapshots which will be used as turbulent field realizations for the propagation of particles. Each snapshot is separated by several cascading timescales at large scales (i.e. ~10 L_{inj}/V, where V is the perturbation velocity). The MHD simulations have been performed at the CINES center on the Jade and Occigen supercalculators. The typical duration of a MHD run at 256^{3}, 512^{3}, and 1024^{3} are 5000, 20 000, and 90 000 h.cpu, respectively.
MHD simulations used in this work.
Before presenting our calculations in detail we want to clarify an issue, barely explicited, concerning the definitions of the Alfvénic Mach number in the literature. In this work we use M_{a} = V/V_{a} with V_{a} defined in the total magnetic field (see Sect. 1) and V is the rms turbulent fluid velocity. Instead, several other works use M_{a0} = V/V_{a0}>M_{a} with V_{a0} defined in the background magnetic field. The turbulence is usually forced in the velocity space so there is additional step to link M_{a} (or M_{a0}) to the level of magnetic fluctuations. Here we consider whereas in many other works η_{0} = δB/B_{0} is used. Table 2 gives the correspondence between η and η_{0}. The advantage of using η is that this quantity is bounded between 0 and 1. In our work both η and M_{a} are averaged over three MHD snapshots and are written ⟨ η ⟩ ≤ ⟨ M_{a} ⟩. In Table 1 the value of ⟨ η ⟩ for the different jobs is given in the fifth column. Comparing these values with those dispalyed in Table 2 we see that our simulations are in the sub to transAlfvénic regime with respect to the Alfvén velocity taken in the background magnetic field.
Correspondence between the parameters η and η_{0}.
2.2. Kinetic module
2.2.1. General description
We solved the Lorentz equation for each particle of charge q and mass m. The particle has a momentum p = γmv and a normalized velocity β = v/c and propagates in an electromagnetic field δE (no mean electric field), B_{T} = δB + B_{0}: (4)The electromagnetic fluctuating components δE and δB are provided by the MHD code. In this work, only protons are considered, hence m = m_{p} and q = + e. Each particle is injected with a Lorentz factor γ_{0}. This defines the particle’s Larmor radius r_{L0} = γ_{0}mc^{2}/eB_{0}. We also define the particle’s synchrotron pulsation Ω_{0} = r_{L0}/c. Equation (4) is normalized with respect to this initial value as (see Beresnyak et al. 2011, hereafter BYL11): We assume that the largescale magnetic field B_{0} is aligned along e_{z}. In the previous equations we have the following notations: , , . In practice, at least in the interstellar medium, the effect of the electric field over highenergy (multi TeV) CRs can be neglected and remains equal to 1.
When particles are submitted to radiative losses (or reacceleration) it is necessary to add another equation for . For instance, relativistic protons cooling under the effect of pion production with a cross section σ_{pp} verify
where t_{pp} ≃ (σ_{pp}n_{H}c)^{1} for a medium of hydrogen density n_{H}. We note that hereafter, as the particle energy is fixed, we denote r_{L0} and Ω_{0} to r_{L} and Ω.
2.2.2. Integration schemes
We tested several integration schemes for Eqs. (5) and (6): a leapfrog scheme (by default implemented in RAMSES to treat the propagation of test particles in a gravitational field), a RungeKutta method of 5th order, and a BulirschStoer method (Press et al. 1992). We tested the different integration methods and found differences at the level of a few percent over the particle mfps. This small difference can be explained by the construction of the code: RAMSES adapts the time step in the PIC module in a way that a particle cannot cross more than one grid cell within one time step (see Teyssier 2002, Sect. 2.4). We used the leapfrog secondorder integration method to save computational resources. The essential effect now is the magnetic field calculation at the particle position (see next section). The typical duration of one run to complete the propagation of 10^{6} particles depends on the particle’s rigidity and the Alfvénic Mach number. It varies from 10^{6} to a few 10^{7} time steps, in other words from ten hours to a few tens of hours of computation on the Jade and Occigen supercomputers at the CINES facility.
2.2.3. Electromagnetic field interpolation
The magnetic field is interpolated at the position of the particle from its values derived by the MHD code at the grid points. For this, we have adopted a volumeaveraged interpolation of the field strength. The volume can include the next eight neighbor points, this is the firstorder cloudincell (CIC) interpolation implemented by default in RAMSES or the next 64 neighbor points using a thirdorder piecewise cubic spline (PCS) interpolation (Haugbølle et al. 2013). The results in this work were obtained using the PCS interpolation method.
3. Direct numerical simulations of cosmic ray transport
Direct numerical simulations of CR transport in magnetic turbulence is a young but already extended field of research. In Sect. 3.1 we summarize the results of two related works proposed by BYL11 and Xu & Yan (2013, hereafter XY13), both developed in the context of CR propagation in the interstellar medium. In Sect. 3.3 we present our results describing CR mfp dependences with respect to Alfvénic Mach numbers while in Sect. 3.5 we present our results describing CR mfp dependences with respect to CR rigidities. In Sect. 3.4 we discuss how forcing geometries affect the CR transport. Finally Sect. 3.6 presents the variations of the ratio of perpendicular to parallel mfps with the particle rigidity at different turbulence levels.
3.1. Summary of previous works
Beresnyak et al. (2011) have investigated the propagation of CRs in incompressible forced MHD turbulence. In the incompressible limit, the effects of the fastmagnetosonic modes are absent, the transport is controlled by the shearAlfvén and the pseudoAlfvén modes (the incompressible limit of slowmagnetosonic modes). Both types of modes have been found to follow a GS scaling (Maron & Goldreich 2001). BYL11 solve Eqs. (5) and (6) neglecting the effect of electric field fluctuations and use different statistically independent (static) magnetic configurations. No backreaction of CRs over the turbulence is considered in this work.
Beresnyak et al. (2011) have conducted 3D simulations at a resolution of 768^{3} in both balanced and imbalanced turbulence cases. Imbalanced turbulence is expected, for instance close to CR sources where a gradient of CR can drive perturbations in a preferred direction with respect to the mean magnetic field or in the solar wind. The turbulence is injected at large scales using a stochastic solenoidal forcing (see the definition in Sect. 2.1.1). The authors have carried out simulations in a regime of weak turbulence (or strong background field) with an Alfvénic Mach number M_{a} = 0.1 and in a regime of strong turbulence (or weak background field) with M_{a} = 1. We note that BYL11 calculate the Alfvén velocity with respect to the background magnetic field B_{0}, hence M_{a} has to be interpreted as M_{a0} = δB/B_{0} in our work. The injection scale of the turbulence in their simulation is (in box length units) L_{inj} ≃ 0.2L. We note that in the subAlfvénic case, BYL11 use an elongated box with a main axis parallel to B_{0}. Then, in that case L_{inj} corresponds to the perpendicular turbulent outerscale and the parallel outerscale of the turbulence in the parallel direction is 10 L_{inj}.
Beresnyak et al. (2011) have computed the variation of spatial diffusion coefficients with respect to the direction of the background magnetic field B_{0}. Balanced and imbalanced turbulence cases do not show strong differences and result in diffusion coefficients varying by a factor of less than a few. In the case of balanced transAlfvénic turbulence, BYL11 find a parallel diffusion coefficient D_{∥} independent of the particle’s rigidity at low rigidity ρ = r_{L}/L ≤ 0.02^{2} then scaling as ρ^{1/3} in the range 0.02 <ρ< 0.1 and scaling as ρ^{2} beyond (see their Fig. 6). Considering perpendicular diffusion, the authors claim to have a coefficient independent of the particle’s rigidity in both subAlfvénic and transAlfvénic cases. In fact, considering their Figs. 7 and 8, their solutions show stronger dependence on rigidity. Especially, D_{⊥} ∝ ρ^{1/3} fits the curve corresponding to the balanced turbulence case in the transAlfvénic case better. We find the same trend in the subAlfvénic case for ρ ≤ 0.05. No explicit effects about the turbulence level has been reported, but a gain of a factor of ~3–4 of difference for D_{⊥} can be deduced between sub and transAlfvénic cases.
Xu & Yan (2013) have investigated the propagation of CRs in compressible MHD turbulence and so the simulations retain the effects of shearAlfvén modes and both magnetosonic modes. Here again the magnetostatic and testparticle limits have been retained.
Xu & Yan (2013) conducted 3D simulations at a resolution of 512^{3}. The turbulence is injected at large scales using a stochastic 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_{a} = 0.19 to M_{a} = 0.73. XY13 have calculated the Alfvén velocity with respect to the total magnetic field B_{T} hence their definition of M_{a} is similar to our. The injection scale of the turbulence in their simulation is (in box length units) L_{inj} = 0.4 L.
Xu & Yan (2013) computed spatial diffusion coefficients with respect to the direction of the background magnetic field B_{0}. The authors conducted a detailed study of the effect of the turbulence level over the particle transport and have calculated CR parallel and perpendicular mfps dependences with respect to M_{a}. The results, however have been presented only at one particle rigidity ρ = 0.01 (so in box length units). Parallel mfps have been found to be >L_{inj} scaling roughly as . Perpendicular mfps have been found to be consistent with a scaling expected theoretically by Yan & Lazarian (2008) in the case the parallel mfps verify λ_{∥}>L_{inj} (see details in Sect. 3.3.3).
We present our results below. We provide an extensive investigation of CR parallel and perpendicular mfps dependences with respect to M_{a} and ρ that extend the two above works significantly. We find powerlaw variations for λ_{∥} and λ_{⊥}. The powerlaw index is denoted α hereafter.
3.2. CR mean free paths: calculation procedure
Cosmic ray mfps are reconstructed using an averaging procedure over a large number of individual particle trajectories. We record the position x,y,z in the 3D simulation cube for each particle: z(t) determines the position of the particle with respect to the background magnetic field B_{0}, x(t),y(t) determines the position of the particle in the directions perpendicular to the background magnetic field B_{0}.
To proceed to mfp calculations, we have selected a set of N_{r} = 3 magnetic field realizations separated by at least two largescale cascade times in order to ensure a statistical independence between two successive realizations. For each realization an ensemble of N_{p} particles are propagated until a convergence is obtained (see next). The parallel spatial diffusion coefficient is calculated using the following formula: (7)The perpendicular spatial diffusion coefficient is calculated using the following formula: (8)These two coefficients are calculated until they both converge to a plateau (see an example in Fig. 1 for a particle at ρ = r_{L}/L = 0.037). We note that at low M_{a} the ballistic regime is present up to t ~ 10^{3}Ω^{1}. Once the diffusion coefficient κ_{∥ , ⊥} = limt → + ∞κ_{∥ , ⊥}(t) is found the corresponding mfp is obtained using the relation λ_{∥ , ⊥} = 3κ_{∥ , ⊥}/v (the particle velocity is identified to the light velocity (v = c) in the cases considered here). For the results presented below a typical number of N_{p} = 10^{6} particles is used.
Fig. 1 Timedependent evolution of λ_{∥}(t) (continuous lines) and λ_{⊥}(t) (dashed lines) for a normalized particle Larmor radius ρ = r_{L}/L = 0.037 at three different Alfvénic Mach numbers with χ = 1 (jobs 9J09c1.0, 9J11c1.0 and 9J12c1.0). The mfp corresponds to the plateau reached at large timescales. Time is in units of particle gyroperiod τ = Ω^{1}. 

Open with DEXTER 
3.3. CR mean free paths: Alfvénic Mach number dependences
We first present the dependence of λ_{∥} and λ_{⊥} with respect to the turbulence level M_{a}. In the following figures, the selected ρ are 0.01,0.037,0.047,0.058,0.074,0.097,0.117,0.147, and 0.186 respectively. We note that all rigidities except for ρ = 0.01 can be associated with a scale in the inertial range of the turbulence as can be seen in Fig. 2.
Another important aspect to keep in mind while comparing these results with the transport in the interstellar medium (ISM) or interplanetary medium (IPM) is that our solutions are restricted to highenergy particles. To fix the ideas we consider typical magnetic fields of 5 μG in the ISM and 50 μG in the solar wind. Hence, considering turbulence injection scales expressed in parsec and astronomical units, we have and in the ISM and IPM, respectively. Here E_{PeV} and E_{GeV} are the particle energies expressed in PeV and GeV units.
Fig. 2 Stationary kinetic energy turbulent spectra corresponding to jobs 8J01c0.5, 9J06c0.5, and 10J01c0.5 in Table 1. Also displayed: the normalized wavenumbers kL corresponding to different normalized Larmor radii ρ = r_{L}/L. 

Open with DEXTER 
3.3.1. Parallel mean free paths
In Fig. 3, we show the dependence of the CR parallel mfps with respect to the Alfvénic Mach number (the upper figure displays the results obtained for χ = 1 (SoF), the lower figure displays the results obtained for χ = 0 (CoF)). It can be seen that parallel mfps decrease with M_{a} at all particle rigidities. The SoF solutions show faster variations with M_{a} than the CoF solutions. We find typical ratios λ_{∥}(χ = 1) /λ_{∥}(χ = 0) varying between ~100 at low M_{a} to ~1 as M_{a} → 1. We find that λ_{⊥}(χ = 1) /λ_{⊥}(χ = 0) varies from 4–5 for M_{a}< 0.6 to ~1 as M_{a} → 1. In Fig. 3 (the SoF case) it can also be seen that the parallel mfp decreases with an increasing rigidity, which is an unsual behavior. We come back to this issue in Sect. 3.5.4.
We extracted a subset of solutions at two different Larmor radius ρ = 0.01 and ρ = 0.097. The former allows us to make some comparisons with XY13, while the latter corresponds to a particle Larmor radius in resonance with a wave number k_{0.097} deep in the inertial range of the turbulence. It corresponds to a scale reasonably far from the injection zone (typically k_{0.097}L ~ 10,k_{0.097}L_{inj} ~ 2), as in our case L_{inj} ~ 0.2L. Particles with rigidities in the inertial range are interesting for investigating the effect of resonant pitchangle scattering over particle parallel mfps. The difficulty is that resonant interactions occurring at kr_{L} ~ 1 are important in the calculation of λ_{∥}, and nonresonant interaction produced by the transittime damping (TTD) mechanism due to fluctuations at scales larger than r_{L} also contribute (Lee & Voelk 1975). The error bars correspond to standard deviations produced by the solutions from the different snapshots.
In the case of incompressible (solenoidal) forcing χ = 1 we find the following solutions: (9)In the case of compressible forcing χ = 0 we find the following solutions: (10)The error on the index is obtained considering the extreme slopes that fit all data points. The mean index values and their errors depend on the errors produced by each data point. We have compared our results using one, three and six snapshots for the job 9J06c0.5. We found that mfps in each of the last two cases are contained within the error of the simulation at one snapshot. The errors obtained using three and six snapshots are similar to (even if slightly reduced in the case with six snapshots), but are reduced with respect to the case with one snapshot. In order to optimize the averaging procedure we keep the number of snapshots N_{r} = 3 in all our calculations.
Fig. 3 Upper figure: parallel mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: parallel mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a compressible forcing (χ = 0). and dependences are shown in the insert. 

Open with DEXTER 
3.3.2. Perpendicular mean free paths
The perpendicular mfps are controlled by two efects: a first contribution is associated with the different processes entering in the parallel transport (resonant and nonresonant pitchangle scattering) and a second contribution is associated with the wandering of magnetic field lines (Jokipii 1971). Figure 4 shows the CR perpendicular mfps as a function of the Alfvénic Mach number.
In the case of solenoidal forcing χ = 1 we find the following solutions: (11)In the case of compressible forcing χ = 0 we find the following solutions: (12)
Fig. 4 Upper figure: perpendicular mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: perpendicular mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a compressible forcing (χ = 0). 

Open with DEXTER 
3.3.3. Discussion
Comparison with previous works: we first compare our results in the SoF case (χ = 1) with the publication of XY13 where the mfp was calculated at ρ = 0.01 only. Comparing our Fig. 3 (upper figure, red points) with Fig. 5 in XY13 we find some discrepancies especially in the regime M_{a}< 0.7 where we find λ_{∥}/L> 100, whereas XY13 has typical mfps λ_{∥}/L< 10. At higher Alfvénic Mach numbers the two results are compatible. This effect is due to the much larger index α ~ −7.7 obtained in our simulations (see Eq. (9)). Our CoF results at ρ = 0.01 are in good agreements with XY13 solutions.
If we now compare the perpendicular mfp results (the red points in our upper Fig. 4 and Fig. 6 in XY13) we find reasonably good agreement: our index is α = 2.94 ± 0.69 and is compatible with the XY13 index α = 4.21 ± 0.75. The normalization factor of the mfps is in agreement since we normalize λ_{⊥} to the box length L, whereas XY13 normalized λ_{⊥} to the injection scale of the turbulence L_{inj}. If the wandering of magnetic field lines dominates the CR perpendicular transport it is not unrealistic to obtain similar results for λ_{⊥} in both studies.
At this stage it is difficult to explain the reason for the discrepancy found between the two works for λ_{∥} at low and moderate M_{a}. What can be said is that even if the forcing is incompressible in both cases, the forcing method is not similar in the two MHD codes. Also, the structure of the two MHD codes and the methods used in the codes to derive fluid solutions are probably different. Another critical aspect is that XY13 derived their solutions at a normalized rigidity ρ = 0.01. Particles at this rigidity have a Larmor radius that is in resonance with perturbations that are in the dissipation range of the turbulence for simulations at a resolution of 512^{3} (as can be clearly seen in Fig. 2). Waveparticle resonance is a critical process that controls particle scattering along the magnetic field lines. It can be seen especially in the SoF case (and to a lesser extent in the CoF case) as all rigidities except for ρ = 0.01 (and also ρ = 0.037 at low M_{a}) have a uniform dependence on M_{a}. This is the main reason why with regard to the problem of propagation of CR in MHD turbulence – we suggest that the results (here and in XY13) derived at ρ = 0.01 have to be interpreted with some care at least concerning λ_{∥}. This is also why hereafter, unless otherwise specified, we concentrate our analysis on Larmor radii that correspond to scales in the inertial range of the turbulent spectrum (that is to say all rigidities but ρ = 0.01 in Fig. 2).
Theoretical frameworks: QLT predicts because in this approximation the scattering frequency is ∝ (see e.g. Sun 2011)^{3}. Casse et al. (2002) in the case of isotropic turbulence have obtained in the limit δB ≪ B_{0}. They have found this scaling to be valid in the regime η → 1, i.e. beyond the QLT validity domain (see their Appendix A)^{4} . A deviation from the QLT result is also expected in the case of a slabtype turbulence with a steep inertial spectrum index^{5} a result obtained using a second order quasilinear theory (Shalchi et al. 2009).
For perpendicular mfps, QLT is obtained in the limit of vanishing scattering (or infinite mfp) along the mean magnetic field (see Jokipii & Parker 1968; Sun 2011). In that regime, field line wandering controls the perpendicular CR transport and leads to . However, recent theoretical predictions have been proposed by Yan & Lazarian (2008; see the discussion in their Sect. 5)^{6} in the framework of the compressible MHD turbulence model. SubAlfvénic turbulence is characterized by two regimes (Lazarian 2006): between L_{inj} and the turbulence is weak, the magnetic power spectrum develops in the perpendicular direction with respect to the mean magnetic field and scales as ; beyond ℓ_{tr} up to scales such that kL_{inj} ~ 40–50, which fall in the dissipative range, GS scaling applies and the turbulence is strong. The CR perpendicular mfp is mostly controlled by the wandering of magnetic field lines but modulated by the effect of scattering: if λ_{∥}/L_{inj}> 1, then perturbations uncorrelated at scales larger than ℓ_{tr} produce a perpendicular mfp . If conversely λ_{∥}/L_{inj}< 1, then diffusion along the field lines lead to . Finally we note that Casse et al. (2002) in the case of isotropic Kolmogorov turbulence numerically obtained λ_{⊥} ∝ η^{4.6}λ_{∥}, i.e. from the above discussion λ_{⊥} ∝ η^{2.6}.
Lazarian & Yan (2014) explored the effect of chaotic magnetic field behavior over CR transport. The authors confirmed that the CR perpendicular transport is diffusive in the regime λ_{∥}/L_{inj}> 1 for curvilinear distances s along the magnetic field lines larger than a few tens of L_{inj}. However, they described the existence of a superdiffusive perpendicular transport regime produced by the strong turbulence regime at intermediary scales s/L_{inj} ~ 1–10 (see Eq. (10) in Lazarian & Yan 2014).
Parallel mean free path: from the upper part of Fig. 3 it can be seen that our SoF results are incompatible with the predictions of QLT as α< −2. Even if we remove the two highest M_{a} points where QLT is questionable, we still find α< −2. We will consider this aspect again in Sect. 3.5.4 while discussing the magnetic power spectra in the SoF limit. In the CoF case, we obtain parallel mfps with indices close to –2 and this even including the highest M_{a} points. This means that our results are consistent with QLT predictions in the low M_{a} regime and extend it in the high M_{a} regime. Our CoF results are consistent with the theoretical predictions of Casse et al. (2002).
Comparing the upper and the lower parts of Fig. 3, it is clear that SoF (χ = 1 ) simulations produce a turbulent spectrum that leads to a less efficient CR scattering (see also Sect. 3.5.4). The SoF favors the production of incompressible Alfvén modes with respect to compressible (in particular fastmagnetosonic) modes even if both types of modes are driven at large scales. At smaller scales, the production of fast modes from Alfvén modes has been found to be weak at low levels of turbulence (Cho & Lazarian 2003). One can expect that the SoF solutions at low M_{a} provide a good hint of the effect of Alfvén waves over the CR transport. Alfvén waves at the lowest perturbation order do not contribute to TTD but only to Landausynchrotron resonance (also known as gyroresonance), so parallel mfps in the SoF case likely result from this process. Recent transport models (see e.g. Chandran 2000) concluded that Alfvénic turbulence does not produce an efficient CR scattering because of the strong anisotropy the perturbations develop towards small scales. At these scales, the resonant condition implies that within a Larmor radius a particle interacts with several uncorrelated perturbations resulting in an inefficient pitchangle scattering. As reported above, we find typical ratios λ_{∥}(χ = 1) /λ_{∥}(χ = 0) ~ 10−100 at low M_{a} depending on the particle rigidity. This partly supports the conclusion of the above analysis, but in the mean time the ratio is not as small as predicted by this model. Yan & Lazarian (2002) have improved Chandran’s model by considering an Alfvénic magnetic tensor derived from numerical calculations of Cho & Lazarian (2003) and found a parallel mfp larger by ~4 orders of magnitude above ~1 GeV, which is closer to our results. It is also possible that a fraction of the energy of the driven velocity field is transferred into fast modes and/or that the gyroresonance effect in Afvénic turbulence is stronger than expected. For ρ< 0.03, the particles are scattered mainly because of the mirroring effect produced by the largest perturbed wavelengths. The differences between SoF and CoF solutions can be used to qualitatively estimate the impact of TTD over the CR transport: in a CoF geometry, the fraction of fastmagnetosonic waves is larger and the TTD contributes more to the scattering of particles. At higher rigidities, i.e. those interacting with turbulent perturbations in the inertial range, we can advance two arguments. First, considering only SoF solutions we see that the parallel mfp is reduced by a factor of ~10 (with respect to the low rigidity regime); this shows that gyroresonance by Alfvén waves can contribute to the pitchangle scattering substantially. Second, if the SoFtoCoF mfp ratio is reduced by a factor of ~10, then the gyroresonance process by Alfvén and fast waves contributes substantially. At higher Alfvénic Mach numbers the differences between the two forcing solutions vanish (the ratio λ_{∥}(χ = 1) /λ_{∥}(χ = 0) → 1 whatever the rigidity regime). It is difficult at this stage to draw conclusion from this effect. It could be associated with a stronger mode conversion between Alfvén and fastmagnetosonic waves at scales smaller than the injection scales.
However, the above analysis is still limited and a quantitative analysis would require to isolate the effects of the different modes that compose the turbulent spectrum in the different forcing limits. One possibility would be to use the mode decomposition procedure proposed by Cho & Lazarian (2003) and improved by Kowal & Lazarian (2010), which should be valid in the limit M_{a}< 1 (see also the cases discussed in Wisniewski et al. 2012).
Perpendicular mean free paths: considering now Fig. 4, we can first compare our results with the solutions obtained by Yan & Lazarian (2008). In the SoF case, λ_{∥}/L_{inj} varies in the range 1–10^{3} (see the upper part of Fig. 3). At ρ = 0.01 we have λ_{∥}/L_{inj}> 100 and . So the result is not compatible with a scaling, but is still close to it. At such low rigidities and because λ_{∥}/L_{inj} ≫ 1, particles are expected to move along a field line over several coherence lengths with a curvilinear abscissa s ~ vt and their diffusion perpendicular mean magnetic field is determined by the divergence of magnetic field lines. Again, because waveparticle interactions occur in the turbulence dissipation range at these rigidities, it is not possible to fully trust that the derived parallel mfp is the one to be expected if the turbulence were fully developed at these scales. Including all the missing modes would likely produce a smaller parallel mfp. In turn, we cannot discard the possibility that particle scattering could modify the result obtained for λ_{⊥} substantially even if we expect the effect to be small at low M_{a}. We do not find any clear trend of superdiffusive CR transport which could produce a larger final perpendicular mfps with respect to the diffusive regime. Superdiffusion is produced by the strongAlfvénic turbulence (Lazarian & Yan 2014) in an interval which widens as M_{a} increases whereas perturbed scales associated with weakAlfvénic turbulence which produces magnetic field line diffusion are restricted to scales ~L_{inj} only. We think that we need a higher grid resolution at low M_{a} in order to include scales where GS spectrum develops in order to capture a possible effect of superdiffusion. At ρ = 0.097, we find , so we have a slower Alfvénic Mach dependence at this rigidity. Again, the result is not compatible with a dependence. At high rigidities there is a competing effect between magnetic field line wandering and particle magnetic scattering. We have λ_{∥}/L_{inj} ~ 1, and not >1, so waveparticle interactions have stronger effect over the perpendicular transport. This may explains why our solutions are not completely compatible with . It is again difficult to isolate an effect of superdiffusion in that rigidity regime. At high M_{a}, we find a perpendicular mean free path about two times smaller than the prediction in relative agreement with the limit found by Yan & Lazarian (2008) and in agreement with the results presented in XY13. In the CoF case, our solution is (lower part of Fig. 4). As λ_{∥} is reduced with respect to the SoF case we also tested the relation . We find that our result is also not completely compatible with the latter relation (although not so far because the extreme values of α are 1.79 and 1.85), respectively, but it is strongly incompatible with . Here again, the parallel mfps cover a range (~ 1−50)L_{inj}. This may explain why we do not find a clear trend. The CoF cases are characterized by a stronger effect of particle scattering produced by fastmagnetosonic waves.
In both forcing geometries our solutions bracket the QLT solution (from above for the SoF solutions and from below for the CoF solutions). Our results are not completely compatible with such a scaling possibly because the last point at M_{a} = 0.89 is a bit under (above) the extrapolation of in the CoF case (SoF case). At high M_{a} indeed the use of QLT is questionable.
Fig. 5 Upper figure: parallel mean free path dependence with M_{a} at ρ = 0.097 for χ = 0,0.5,1. Lower figure: perpendicular mean free path dependence with M_{a} at ρ = 0.097 for χ = 0,0.5,1. 

Open with DEXTER 
Fig. 6 Upper figure: parallel mfps at three different resolution for jobs 8J02c0.5, 9J06c0.5, and 10J01c0.5. Lower figure: the same but for perpendicular mfps. 

Open with DEXTER 
Summary: our solutions are globally consistent with the results obtained by XY13 to the notable exception of λ_{∥} at M_{a}< 0.7. This effect is possibly related to differences in forcing procedures and methods used to solve the MHD equations in the codes. If we accept that SoF (CoF) drives at large scales Alfvén (fastmagnetosonic) waves preferentially, our results can give some qualitative hints about the respective effect of gyroresonance and TTD processes. At low rigidity, i.e. for particles in resonance with modes in dissipation range of the turbulence, 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 resolution of the simulations. With respect to the low rigidity findings, at higher rigidities the mfp is considerably reduced in the SoF case, which can be due to the effect of gyroresonance by the Alfvénic turbulence. Also, the ratio of the parallel SoF and CoF mfps decreases, which can be due to the effect of gyroresonance of Alfvén and fastmagnetosonic waves over the CR transport. We find that predicted by the QLT is compatible with CoF solutions, but not with with SoF solutions. Our results concerning λ_{⊥} are less affirmative than the ones advanced by XY13 and do not show a clear trend at the moment. The results have not been found to be fully compatible with the theoretical predictions proposed in Yan & Lazarian (2008). One explanation is that for each specific χ value the parallel mfp is not entirely in the regime λ_{∥}/L_{inj}< 1 or λ_{∥}/L_{inj}> 1. Finally, our solutions are between the two scalings: and (the latter expected from QLT). We did not find a clear trend in our results that could be associated with an effect of superdiffusion especially at low rigidities and low Alfvénic Mach numbers likely because GS turbulence is restricted to reduced scale lengths. More simulations with larger dynamics are needed in order to choose definitively from among the effects of particle scattering and any effect of magnetic field line transport and also to decide among the different theoretical models. In particular, results at M_{a} ~ 0.1 would greatly help to probe the QLT predictions further. Unfortunately, such low Afvénic Mach numbers require very long integration timescales to reach a diffusive plateau. We come back to this aspect in Sect. 3.5.4.
3.4. CR mean free paths: forcing effects
Figure 5 presents the dependence of parallel and perpendicular mfps on M_{a} at ρ = 0.097 for χ = 0,0.5, and 1.
It can be seen that in both cases (parallel and perpendicular mfps), especially for M_{a} ≤ 0.5, the MoF solutions at χ = 0.5 are closer to the SoF solutions. This is supported by the value of r_{c} and the fraction of compressible modes (see Sect. 2.1.1): in the MoF case r_{c} = 1/3 is closer to r_{c} = 0 obtained for χ = 1 than r_{c} = 1 obtained for χ = 0. This fits with the arguments advanced above that Alfvén waves are preferentially produced with respect to fast magnetosonic waves in the SoF and MoF geometries and that a tiny fraction of the Alfvén energy is transferred to fast waves at low turbulent levels (Cho & Lazarian 2003). At higher Alfvénic numbers the effect of the forcing geometry is less stringent possibly because of a more efficient mode conversion. There, we find that parallel and perpendicular mfps are independent of χ.
3.5. CR mean free paths: Rigidity dependences
3.5.1. Resolution tests
Before discussing any specific rigidity dependence we first test our calculation procedure for ρ> 0.037 at different grid resolution. Figure 6 presents the parallel and perpendicular mfp with respect to ρ at X = 8, 9, and 10. All the simulations were performed at χ = 0.5 and jobs 8J02c0.5, 9J06c0.5, and 10J01c0.5 were used.
It can be seen that the results are compatible each other whatever the value of X. We also see a trend with a decreasing parallel mfp as the resolution increases, especially at rigidities ρ< 0.1. An effect that is expected since as X increases the inertial turbulence zone is shifted towards smaller scales.
3.5.2. Parallel mean free paths
Figure 7 presents CR parallel mfps dependences with respect to the particle normalized Larmor radius ρ = r_{L}/L.
Fig. 7 Upper figure: parallel mfps versus particle normalized Larmor radii ρ at different M_{a} in the SoF case (χ = 1). Insets: slopes obtained for different indices. Lower figure: parallel mfp versus particle normalized Larmor radius ρ at different M_{a} in the CoF case (χ = 0). 

Open with DEXTER 
Parallel mean free path indices and errors for SoF and CoF geometries.
In Table 3 we present the fit indices and the errors corresponding to linear fits of the parallel mfp for 0.01 <ρ< 0.1 and 0.03 <ρ< 0.1 in both SoF and CoF cases respectively.
3.5.3. Perpendicular mean free paths
Figure 8 presents the CR perpendicular mfp with respect to the particle normalized Larmor radius ρ = r_{L}/L.
Fig. 8 Upper figure: perpendicular mfp versus particle normalized Larmor radius ρ at different M_{a} values in the SoF case (χ = 1). Lower figure: perpendicular mfp versus particle normalized Larmor radius ρ at different M_{a} values in the CoF (χ = 0). 

Open with DEXTER 
Perpendicular mean free path indices and errors for both forcing geometries.
In Table 4 we present the fit indices and the errors corresponding to linear fits of the perpendicular mfp for 0.01 <ρ< 0.1 and 0.03 <ρ< 0.1 in the SoF and CoF cases respectively.
3.5.4. Discussion
Some obvious points can be highlighted. The SoF and CoF solutions for λ_{∥} show very different rigidity variations. In the SoF case we see a clear trend that there is a spectral softening as M_{a} decreases with some solutions showing inverted spectra (negative α). In the CoF case no such effect is present and all the solutions have a spectral index largely independent of the Alfvénic Mach number. Solutions for λ_{⊥} show similar spectra in the SoF and CoF cases, even if the SoF solutions have a stronger evolution with M_{a} (see Sect. 3.3.3).
Comparison with previous works: we first compare our results with the solutions obtained by BYL11. As reported above, BYL11 have derived their results in the incompressible MHD limit; that is in the limit β_{p} = P_{g}/P_{m} ~ (M_{a}/M_{s})^{2} → ∞^{7}. As can be seen from Table 1, most of our simulations are performed in an equipartition regime with β_{p} ~ 1. BYL11 also performed simulations at M_{a} = 1 and M_{a} = 0.1, but presented their results for λ_{∥} at M_{a} = 1 only (see Sect. 3.1). Hence, the more appropriate simulation to be compared with these results is the job 9J12c1.0. In that case, we find an index compatible with α = 1/3 and α = 0 (thus retaining or not the point at ρ = 0.01) hence compatible with the BYL11 results (see Table 3). It should be noted here that the results at rigidities ρ< 0.02 presented in BYL11 correspond to particles propagating in the dissipation regime of the MHD turbulence even if the resolution of the MHD simulations performed by the authors is higher than ours (768^{3} instead of 512^{3}). Again, we think these results should be considered with some caution especially for λ_{∥}. For λ_{⊥}, λ_{⊥} ∝ ρ^{1/3} at M_{a} = 0.89 can provide a reasonable fit of the points with ρ< 0.1 in the SoF case but still removing the point at ρ = 0.01. Hence it is also compatible with BYL11 results. In order to push ahead the comparison with BYL11, we also performed a series of simulations at level 8 at high β_{p} ~ 100 (job 8J02c0.1) and M_{a} = 1. Here the index is found to be compatible (although a bit larger) with α = 1/3, but not with α = 0.
Theoretical frameworks: QLT predicts λ_{∥} ∝ ρ^{2−s}, where s is the 1D index of the magnetic turbulent spectrum. In the case of a Kraichnan or a Kolmogorov spectrum with s = 3/2 and s = 5/3, QLT predicts α = 1/2 and 1/3, respectively. However, such indices are not restricted to QLT solutions. For instance, solutions of CR transport in isotropic turbulence with a Kolmogorov spectrum have an index α = 1/3 irrespective of the turbulence level (Casse et al. 2002), a result that cannot be anticipated using the predictions of the QLT. One important question is wether this scaling is also recovered in a GoldreichSridhar (GS) type turbulence (Marcowith et al. 2006). Indeed the GS spectrum follows a Kolmogorov scaling in the perpendicular wave number space. Hussein et al. (2015) performed simulations using different synthetic turbulence models and in particular the GS model. Their simulations have been performed in the strong turbulence limit with η_{0} = δB/B_{0} = 1. The authors found α ≃ 1/3–1/2 (the fit over one decade in rigidity does not allow a value between the two scalings) and also found that GS turbulence is almost as efficient as the isotropic Kolmogorov turbulence to scatter CRs. Yan & Lazarian (2002) and Chandran (2000) found several orders of magnitude of differences between the CR mfps produced by isotropic and GS models for M_{a}< 1^{8}. Yan & Lazarian (2008) instead reported on a dominance of fastmagnetosonic modes which follow an isotropic Kraichnan turbulent spectrum. Hence, if this result is correct, one can expect to find α = 1/2. We note that this result has been derived using a nonlinear correction to the particle trajectory and hence it supersedes QLT solutions (see Yan & Lazarian 2008 for details). Finally, Hussein et al. (2015) obtained λ_{⊥} ∝ ρ^{1/2} in the case of GS spectrum at δB = B_{0} (even if the fit is restricted to almost one decade in rigidity). At lower M_{a}, only if λ_{∥}<L_{inj}Yan & Lazarian (2008) find a perpendicular mfp varying with the particle rigidity as .
Parallel mean free paths: we now discuss each specific result in details. In all cases, the indices for the parallel mfp are found in the range (− 1.5,1). In particular, spectra harder than λ_{∥} ∝ ρ are rejected.
The lower part of Fig. 6 presents the parallel mfp in the CoF case. In that case, we never find indices compatible with the Bohm scaling (α = 1) but the index at M_{a} = 0.89 which is marginally compatible. The mfp index in the regime ρ< 0.1 shows some hardening from low to high M_{a}, except for in all cases indices are compatible either with 1/3 or 1/2 (see Table 3). This hardening trend is reinforced if we include the point at ρ = 0.01 in the fit: the points at ρ = 0.01 produce a bias towards smaller values of the indices. However, it is expected that if located in the inertial range of the turbulence, these points should be shifted towards smaller mfps. Overcoming this effect would have required simulations at levels larger than X = 10, which is beyond the computing resources available for this work.
If we now turn to the SoF case (upper part of Fig. 7) at low M_{a}, the results are clearly incompatible with α = 1/2 or even 1/3 (see Table 3), and either when including or removing the points at ρ = 0.01. For M_{a} ≤ 0.5 negative indices are obtained. This result is difficult to anticipate with the actual theoretical transport models. One possibility would be to invoke intermittency effects known to reduce the particle mean free path with respect to homogeneous turbulence (AlouaniBibi & le Roux 2014). However, this possibility is questionable since we did not detect any subdiffusive paths at low rigidities and low M_{a}. Inverted indices can be associated with a propagation of CRs in a damped turbulence (Yan & Lazarian 2008; Shalchi & Büsching 2010). Above M_{a} = 0.5 our results are compatible with α = 1/3 which is the expected result in Kolmogorov isotropic turbulence (Casse et al. 2002).
From the above discussion we find that λ_{∥} ∝ ρ^{1/2} is clearly compatible with the simulations in the CoF case. Hence we find that in the CoF limit our results are consistent with a CR transport controlled by fastmagnetosonic waves as interpreted in the paradigm presented by Yan & Lazarian (2008). This argument is similar to the conclusion reached from the Alfvénic Mach number analysis in Sect. 3.3.3. However, our results are also consistent with α = 1/3 and a definite answer will require more intensive simulations with extended dynamical range in order to select bewteen the two models. In Fig. 9, we plot the magnetic energy density power spectrum at two Alfvénic Mach numbers M_{a} = 0.34 and M_{a} = 0.89. We see that the power spectrum index is almost the same (s ~ 1.5–1.6), which can explain why we find a rigidity index α = 1/2−1/3 largely independent of M_{a}.
Fig. 9 Magnetic energy density power spectra for M_{a} = 0.88 (blue continuous line) and M_{a} = 0.34 (red continuous line) in the CoF case. 

Open with DEXTER 
The simulations obtained in the SoF case show very different trends. One way to get inverted spectra would be to have s ≥ 2, hence the results could still be interpreted in the QLT framework. We plot in Fig. 10 the magnetic energy density power spectrum at M_{a} = 0.34 and M_{a} = 0.89. We checked that the stationarity was reached. If the spectrum at M_{a} = 0.89 has an index s ~ 1.5 the spectrum at M_{a} = 0.34 is very different: it is harder at large scales with s ~ 1 and softer above kL ~ 18 to s ~ 2. We note that scales kL ~ 18 correspond to rigidities ρ ~ 0.05 in Fig. 2 and this softening effect can explain the flat or negative values of α. However, we cannot exclude with our simulations at level 9 that this spectral effect is of numerical origin; because of the finite resolution our simulations may have missed some resonant modes. In order to test this possibility we performed one simulation at M_{a} = 0.34, but at level 10 (see Table 3). At X = 10 we also find an inverted spectrum for the parallel mfp with an index compatible with the solution obtained at level 9, but smaller in absolute value (the effect is more pronounced if we include the point at ρ = 0.01). This can be expected because low rigidity points with ρ = 0.01,0.037 and 0.047 have lower parallel mfp at X = 10. A higher resolution shifts the dissipation scale towards smaller scales and produces more fluctuations in resonance with low rigidity particles (see Sect. 3.5.1). So, we conclude that inverted spectra probably do not result from a numerical artifact. It would be interesting to explore the (M_{a},M_{s}) parameter space to test if this kind of spectrum is a general trend. Finally, this effect can possibly be related with the behavior of weak subAlfvénic turbulence (see Sect. 3.3). Hence, if M_{a} = 0.3, weak subAlfvénic turbulence is present over one magnitude in k beyond the injection scale. Because the energy cascades in the perpendicular direction, in that case, a low CR scattering efficiency and large parallel mfps are expected.
Fig. 10 Magnetic energy density power spectra for M_{a} = 0.88 (blue continuous line) and M_{a} = 0.34 (red continuous line) in the SoF case. 

Open with DEXTER 
When considering parallel mfps above ρ = 0.1, in all cases the spectrum hardens but our simulations do not have enough dynamics to probe the highrigidity regime where λ_{∥} ∝ ρ^{2} is expected. This rigidity dependence is expected in the configuration of CR propagation with Larmor radii larger than the turbulence coherence length (note that L_{inj} ~ 0.2) (see Deligny et al. (2004), Plotnikov et al. (2013)). We also note that a transition between small to large scale turbulence at ρ = 0.1 gives a transition Larmor radius r_{L} ~ 0.5 L_{inj} as L_{inj} = 0.2L.
Perpendicular mean free paths: considering now perpendicular mfps (see Figs. 7), the CoF solutions are not compatible with α = 1/3, but are compatible with α = 1/2. The SoF solutions show spectral indices all compatible with α = 1/2 with or without the points at ρ = 0.01. In fact, if we remove these points the results are also compatible with α = 1/3 (this may be less true at low M_{a}). At low M_{a}, the index is also compatible with the solution obtained at level 9 (see Table 4). At M_{a} = 0.66 and M_{a} = 0.89, we find that s ~ 1.6 fits the data well leading to α = 1/3, but the spectral dynamics is too restricted again to separate solutions with α = 1/2 or 1/3. At M_{a} ~ 0.89, we also find that the mfps in the CoF and SoF cases are similar. The interpretation of this is difficult and require requires distinguishing the impact of each type of mode (see Sect. 3.3.3). It is likely that λ_{⊥} at high M_{a}, whatever the forcing, is controlled by the field line wandering process, which is an issue difficult to test using an analytical approach (see discussions in Casse et al. 2002; Shalchi 2010). We note that the wandering of the field lines is controlled by the shearing effect produced during the interaction of two counterpropagating Alfvén wave packets. The perpendicular mfp above ρ = 0.1 is compatible with a flat dependence on ρ (see Casse et al. 2002; Hussein et al. 2015).
Summary: we again find that CoF solutions are compatible with QLT predictions (at least at low M_{a}) but also with nonlinear models of CR diffusion in fastmagnetosonic turbulence. However, the dynamics in rigidity scales is not sufficient to differentiate between solutions with α = 1/3 or α = 1/2. SoF solutions, however do not fit in this scenario; they are incompatible with QLT solutions at low M_{a} and show spectra with a milder or inverted rigidity dependence that are connected with particular magnetic field spectra possibly connected with the transition between weak and strong subAlfvénic turbulence regimes. Perpendicular mfps show similar rigidity dependences in both SoF and CoF cases. There the wandering of field lines is likely important even if not dominant.
3.6. The λ_{⊥}/λ_{∥} ratio
Figure 11 presents the ratio λ_{⊥}/λ_{∥} with respect to ρ derived at M_{a} = 0.34,0.66 and 0.89 in the SoF and CoF cases. This diagnostic is interesting for probing the importance of field line wandering over the CR transport (see Casse et al. (2002)). To that end it is first instructive to compare the latter ratio with the prediction from classical scattering theory (Forman & Gleeson 1975): (13)The solution given by Eq. (13) is displayed as the continuous line in the two panels in Fig. 11. In both the SoF and CoF cases Eq. (13) is not compatible with our results (maybe with the exception of highrigidity points where the contribution of pitchangle scattering to λ_{⊥} should be the strongest). The difference amounts to the contribution of the field line wandering process over the CR perpendicular transport. If we remove the point at ρ = 0.01 our results are compatible with a flat ratio for ρ< 0.1 in the CoF case and in the SoF but only at M_{a} = 0.89. If χ = 1 and M_{a}< 0.89 a flat dependence is disfavored. Above ρ = 0.1 the dynamics is reduced but the ratio is compatible with ρ^{2} especially in the CoF case. A flat ratio is expected if the field line wandering process takes over particle scattering. It seems that a reduced efficiency of the parallel transport in subAlfvénic turbulence could provide an explanation of this spectral hardening at low CR rigidities. A more quantitative discussion on field line wandering requires a reconstruction of the magnetic field lines in different MHD realizations. This is the subject of a work in preparation.
Fig. 11 Ratios of λ_{⊥}/λ_{∥} for different values of the Alfvénic Mach number. Upper figure: SoF case. Lower figure: CoF case. The continuous curves display the function 1 /^{(}1 + (λ_{∥}/r_{L})^{2}^{)}. 

Open with DEXTER 
4. Summary and conclusions
In this work we have developed a series of kineticMHD simulations of CR transport in magnetized turbulence. We have upgraded the RAMSES code to include a module that permits forcing the velocity space following different geometries and a particleincell module which allows the reconstruction of the trajectory of charged particles in an electromagnetic field. We have investigated the Alfvénic Mach number and rigidity dependence of both parallel and perpendicular mfps by sampling a large number of particles and averaging over several statistically independent field realizations.
The main results of this work are the following:

Forcing effects: The effects of the forcing geometry are the strongest at low Alfvénic Mach numbers where SoF produce parallel mfp that can be two orders of magnitude larger than CoF mfp solutions. A MoF produces results close to the SoF solutions. At high Alfvénic Mach numbers, the impact of the forcing geometry vanishes and all the results converge to a similar solution.

Alfvénic Mach number dependence: In the SoF case, we found faster M_{a} dependences with respect to the results of XY13, but our results are compatible at M_{a}> 0.7. Even if we found a similar trend to XY13 for the dependence of λ_{⊥} on M_{a}, we consider that these results have to be taken with care because these correspond to particle rigidities lying in the dissipation range of the turbulence. The missing resonant interactions are expected to reduce λ_{∥} and possibly modify the result on λ_{⊥}. SoF results are not found to be compatible with QLT predictions. In the CoF case, we found a scaling compatible with and so compatible with QLT predictions but for an Alfvénic Mach number regime beyond the standard validity domain of the QLT. SoF perpendicular mfps are not found to be compatible with the theoretical calculations proposed by YL08, even if our SoF solutions are close to a scaling. This can be explained because the ratio λ_{∥}/L_{inj} is neither ≫ 1 nor ≪ 1 in our solutions. We did not find a clear signature of superdiffusion in the CR perpendicular mfp produced by GStype turbulence possibly because of the limited resolution level. CoF perpendicular mfps have scalings relatively close to the QLT expectation owing to the wandering of magnetic field lines.

Rigidity dependence: CoF gives results that are compatible with λ_{∥} ∝ ρ^{1/2}, consistent with predictions from QLT and the diffusion produced by a wave spectrum following a Kraichnan law. This solution is expected if fastmagnetosonic waves dominate the CR scattering process and is valid even beyond QLT calculations. However, our results are also compatible with λ_{∥} ∝ ρ^{1/3}. Simulations with larger dynamics are necessary in order to differentiate bewteen the two scalings. SoF show very different results, in particular at low Alfvénic Mach numbers where inverted or flat spectra are obtained. This result cannot be explained by QLT predictions. The effect may result from an inefficient CR scattering in weak subAlfvénic turbulence. Perpendicular mfps have rigidity dependences also compatible with ρ^{1/2} for both CoF and SoF cases. CoF and SoF solutions are similar at high M_{a}. At rigidities ρ> 0.1, that is at rigidities larger than the turbulence injection scale, a hardening (a softening) of the parallel (perpendicular) mean path is obtained compatible with the propagation of particles in smallscale turbulence.

Ratio of perpendicular to parallel mfps: A flat ratio characteristic of the effect of magnetic field line wandering at rigidities ρ< 0.1 is found in both CoF and SoF cases at high M_{a}. CoF solutions show the same trend whatever the level of turbulence whereas SoF solutions due to enhanced parallel mfp show a rigidity dependent ratio at low Alfvénic Mach number. Here weak scattering effects observed for λ_{∥} seem to be as important as field line wandering to explain a smaller ratio at low rigidities.
Several issues emerge from this work. The solutions are found to be dependent on the adopted forcing geometry. It appears that CoF simulations are compatible with expectations from QLT or from some nonlinear models, whereas SoF solutions are not. We have proposed some arguments to explain the parallel mfp in the SoF and CoF geometries. They are based on a preferential mode type production and the relative effect of TTD and gyroresonance. However, this is speculation and it seems important at this stage to have a proper test of the effect of each type of MHD mode on the CR scattering process. This approach is important, for instance, to assert the inverted spectra obtained at low Alfvénic Mach numbers in the SoF case. To that end, it would be interesting to filter the different modes out of the magnetic field realizations and to repeat particle propagation runs. This aspect deserves a future investigation. This issue is also connected to the scattering efficiency in an Alfvénic turbulence which has been advanced to be weak because of the anisotropy but has never been tested by direct numerical simulations. It also seems important to extend the calculations of magnetic field transport in the low M_{a} ≤ 0.1 Alfvénic Mach number regime. It would be also interesting to perform highresolution level simulations with X beyond 10 in order to be able to capture and disentangle the effect of strong and weak turbulence regimes over CR perpendicular mean free paths. These calculations however, require dedicated important computational resources that deserve future work. Extension of the present study to the trans and superAlfvénic (M_{a} ≥ 0.9) cases would also be interesting in order to test the different models describing the propagation of particles in isotropic turbulence. Finally, it is also important to extend the simulations at low and high β_{p} regimes as most of the jobs in this work were performed at β_{p} ~ 1. The magnetic field line transport process is important in order to understand the perpendicular mfp of the particles. Dedicated calculations of field line calculations in different forcing geometries are in progress and deserve a forthcoming work.
Finally, it should be emphasized that this study is restricted to the transport of particles in largescale turbulence which is, for instance the case of highenergy CRs (TeV and beyond) in the ISM. The transport of lowenergy particles in the ISM requires adding a component of selfgenerated waves (see the discussion in Sect. 6 of XY13). As the pressure imparted in the lowenergy CR is as important as the gas and magnetic pressure, this would be possible using PICMHD simulations, but only if the source terms associated with the CR currents are correctly implemented in the MHD equations (see Bai et al. 2015). The corresponding upgrading of the PIC and MHD modules is inprogress.
and k are the perturbed wave numbers defined along the local magnetic field B_{L} and the total magnetic field B_{T}, respectively. We also define the global or background magnetic B_{0}, i.e. the large scale background magnetic field. The total magnetic field includes the contribution of the background magnetic field and the perturbed magnetic field δB. The local magnetic field is not unique. Its direction varies with the scale under specification (see the discussion in Cho & Vishniac 2000). ℓ is the coherence length of the turbulence.
Jokipii & Parker (1968) instead predicted and .
Casse et al. (2002) defined η as a quantity which corresponds to η^{2} in our text.
This is the 1D wave number index, denoted s in Sect. 3.5.4. Deviations from QLT are obtained for s> 2.
Yan & Lazarian (2008) used M_{a0} instead of M_{a}.
Aalthough the mfps obtained by Yan & Lazarian (2002) are 4 orders of magnitude above the ones obtained by Chandran (2000).
Acknowledgments
We thank M. Lemoine, P. Hennebelle and H.Yan for fruitful discussions. We thank the anonymous referee for his/her references to the weak turbulence regime and superdiffusive effects. The simulations have been performed using the super calculator facilities available at the CINES (Centre d’Informatique de l’Enseignement Supérieur in Montpellier) center. This work has benefit from the support of the ANR COSMIS project.
References
 AlouaniBibi, F., & le Roux, J. A. 2014, ApJ, 781, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Beresnyak, A., Yan, H., & Lazarian, A. 2011, ApJ, 728, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Candia, J., & Roulet, E. 2004, J. Cosmol. Astropart. Phys., 10, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Casse, F., Lemoine, M., & Pelletier, G. 2002, Phys. Rev. D, 65, 023002 [NASA ADS] [CrossRef] [Google Scholar]
 Chandran, B. D. G. 2000, Phys. Rev. Lett., 85, 4656 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J., & Lazarian, A. 2003, MNRAS, 345, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J., & Vishniac, E. T. 2000, ApJ, 539, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Deligny, O., LetessierSelvon, A., & Parizot, E. 2004, Astropart. Phys., 21, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Farmer, A. J., & Goldreich, P. 2004, ApJ, 604, 671 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., RomanDuval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Federrath, C., Chabrier, G., Schober, J., et al. 2011, Phys. Rev. Lett., 107, 114504 [NASA ADS] [CrossRef] [Google Scholar]
 Forman, M. A., & Gleeson, L. J. 1975, Ap&SS, 32, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Haugbølle, T., Frederiksen, J. T., & Nordlund, A. 2013, Phys. Plasmas, 20, 062904 [NASA ADS] [CrossRef] [Google Scholar]
 Hussein, M., Tautz, R. C., & Shalchi, A. 2015, J. Geophys. Res. (Space Phys.), 120, 4095 [NASA ADS] [CrossRef] [Google Scholar]
 Jokipii, J. R. 1966, ApJ, 146, 480 [NASA ADS] [CrossRef] [Google Scholar]
 Jokipii, J. R. 1971, Rev. Geophys. Space Phys., 9, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Jokipii, J. R., & Parker, E. N. 1968, Phys. Rev. Lett., 21, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Kowal, G., & Lazarian, A. 2010, ApJ, 720, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Laitinen, T., Dalla, S., & Kelly, J. 2012, ApJ, 749, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A. 2006, ApJ, 645, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Lazarian, A., & Yan, H. 2014, ApJ, 784, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. A., & Voelk, H. J. 1975, ApJ, 198, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Marcowith, A., Lemoine, M., & Pelletier, G. 2006, A&A, 453, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maron, J., & Goldreich, P. 2001, ApJ, 554, 1175 [NASA ADS] [CrossRef] [Google Scholar]
 Plotnikov, I., Pelletier, G., & Lemoine, M. 2013, MNRAS, 430, 1280 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing (New York: Cambridge University Press) [Google Scholar]
 Ptuskin, V. S., Moskalenko, I. V., Jones, F. C., Strong, A. W., & Zirakashvili, V. N. 2006, ApJ, 642, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Reville, B., O’Sullivan, S., Duffy, P., & Kirk, J. G. 2008, MNRAS, 386, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin: Springer) [Google Scholar]
 Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shalchi, A. 2009, Nonlinear Cosmic Ray Diffusion Theories, Astrophys. Space Sci. Libr., 362 [Google Scholar]
 Shalchi, A. 2010, ApJ, 720, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A., & Büsching, I. 2010, ApJ, 725, 2110 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A., Koda, T. Å., Tautz, R. C., & Schlickeiser, R. 2009, A&A, 507, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sun, P. 2011, International Cosmic Ray Conference, 10, 240 [NASA ADS] [Google Scholar]
 Tautz, R. C. 2010, Comput. Phys. Comm., 181, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Völk, H. J. 1973, Ap&SS, 25, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Wisniewski, M., Spanier, F., & Kissmann, R. 2012, ApJ, 750, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Wisniewski, M., Kissmann, R., Spanier, F., & Spanier. 2013, J. Plasma Phys., 79, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, S., & Yan, H. 2013, ApJ, 779, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2002, Phys. Rev. Lett., 89, 1102 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2008, ApJ, 673, 942 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2011, ApJ, 731, 35 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Timedependent evolution of λ_{∥}(t) (continuous lines) and λ_{⊥}(t) (dashed lines) for a normalized particle Larmor radius ρ = r_{L}/L = 0.037 at three different Alfvénic Mach numbers with χ = 1 (jobs 9J09c1.0, 9J11c1.0 and 9J12c1.0). The mfp corresponds to the plateau reached at large timescales. Time is in units of particle gyroperiod τ = Ω^{1}. 

Open with DEXTER  
In the text 
Fig. 2 Stationary kinetic energy turbulent spectra corresponding to jobs 8J01c0.5, 9J06c0.5, and 10J01c0.5 in Table 1. Also displayed: the normalized wavenumbers kL corresponding to different normalized Larmor radii ρ = r_{L}/L. 

Open with DEXTER  
In the text 
Fig. 3 Upper figure: parallel mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: parallel mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a compressible forcing (χ = 0). and dependences are shown in the insert. 

Open with DEXTER  
In the text 
Fig. 4 Upper figure: perpendicular mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: perpendicular mfp versus Alfvénic Mach numbers M_{a} at different particle rigidities in the case of a compressible forcing (χ = 0). 

Open with DEXTER  
In the text 
Fig. 5 Upper figure: parallel mean free path dependence with M_{a} at ρ = 0.097 for χ = 0,0.5,1. Lower figure: perpendicular mean free path dependence with M_{a} at ρ = 0.097 for χ = 0,0.5,1. 

Open with DEXTER  
In the text 
Fig. 6 Upper figure: parallel mfps at three different resolution for jobs 8J02c0.5, 9J06c0.5, and 10J01c0.5. Lower figure: the same but for perpendicular mfps. 

Open with DEXTER  
In the text 
Fig. 7 Upper figure: parallel mfps versus particle normalized Larmor radii ρ at different M_{a} in the SoF case (χ = 1). Insets: slopes obtained for different indices. Lower figure: parallel mfp versus particle normalized Larmor radius ρ at different M_{a} in the CoF case (χ = 0). 

Open with DEXTER  
In the text 
Fig. 8 Upper figure: perpendicular mfp versus particle normalized Larmor radius ρ at different M_{a} values in the SoF case (χ = 1). Lower figure: perpendicular mfp versus particle normalized Larmor radius ρ at different M_{a} values in the CoF (χ = 0). 

Open with DEXTER  
In the text 
Fig. 9 Magnetic energy density power spectra for M_{a} = 0.88 (blue continuous line) and M_{a} = 0.34 (red continuous line) in the CoF case. 

Open with DEXTER  
In the text 
Fig. 10 Magnetic energy density power spectra for M_{a} = 0.88 (blue continuous line) and M_{a} = 0.34 (red continuous line) in the SoF case. 

Open with DEXTER  
In the text 
Fig. 11 Ratios of λ_{⊥}/λ_{∥} for different values of the Alfvénic Mach number. Upper figure: SoF case. Lower figure: CoF case. The continuous curves display the function 1 /^{(}1 + (λ_{∥}/r_{L})^{2}^{)}. 

Open with DEXTER  
In the text 