EDP Sciences
Free Access
Issue
A&A
Volume 588, April 2016
Article Number A73
Number of page(s) 14
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201527376
Published online 21 March 2016

© 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. long-wavelength perturbations with scales comparable to the particle’s Larmor radius rL = / Ω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 quasi-linear 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 pitch-angle particle scattering and particle distribution isotropization times. The main drawback of the QLT is that it leads to infinite cosine pitch-angle (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 non-linear solutions which have to be retained in diffusion coefficient calculations, are model dependent. Needless to say that all these studies are performed in the test-particle 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 slow-magnetosonic modes have been found to follow the same GS scaling, but fast-magnetosonic 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 above-mentioned 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 model-dependent) 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 non-relativistic (and also relativistic) shocks where particle-in-cell (PIC) techniques calculate the electromagnetic field solutions from the Maxwell equations. Recent PIC simulations have started self-consistent 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 PIC-MHD 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 sub-Alfvénic turbulence cases. We note that in this work the Afvénic Mach number Ma is defined as the ratio of the rms turbulent fluid velocity V = ⟨ δV21/2 to the Alfvén speed taken with the total magnetic field with and where ρ is fluid mass density. Hence hereafter Ma< 1, but also δB/B0 ≤ 1 so that we can always define a parallel and a perpendicular transport with respect to the background magnetic field B0.

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 set-up 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. Kinetic-MHD 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. Magneto-hydrodynamic 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 Ornstein-Uhlenbeck 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 Nm modes (typically Nm = 32). One advantage of the Ornstein-Uhlenbeck 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 auto-correlation 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 Ornstein-Uhlenbeck process. Hence, χ = 0 corresponds to a pure compressive forcing (hereafter CoF or curl-free) and χ = 1 corresponds to a pure solenoidal forcing (hereafter SoF or divergence-free). 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 (right-hand side terms in Eq. (3)). One gets rc = Fcomp/ (Fsol + Fcomp) = (1 − χ)2/ (1−2χ + 3χ2) (Federrath et al. 2010). Hence rc = 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, large-scale gravitational contraction, supernova explosions, and proto-stellar 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 set-up are P = ρ = 1 unless otherwise specified and the sound speed is . The energy is injected into the large-scale 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 = 2X. We note that hereafter we make the distinction between L, which is the size of simulation cube, and Linj<L, which is the scale of turbulence injection identified with the coherence length of the turbulence.

2.1.2. Simulation set-up

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 = 5123 = 29 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 quasi-stationary 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 Linj/V, where V is the perturbation velocity). The MHD simulations have been performed at the CINES center on the Jade and Occigen super-calculators. The typical duration of a MHD run at 2563, 5123, and 10243 are 5000, 20 000, and 90 000 h.cpu, respectively.

Table 1

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 Ma = V/Va with Va defined in the total magnetic field (see Sect. 1) and V is the rms turbulent fluid velocity. Instead, several other works use Ma0 = V/Va0>Ma with Va0 defined in the background magnetic field. The turbulence is usually forced in the velocity space so there is additional step to link Ma (or Ma0) to the level of magnetic fluctuations. Here we consider whereas in many other works η0 = δB/B0 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 Ma are averaged over three MHD snapshots and are written η ⟩ ≤ ⟨ Ma. 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 trans-Alfvénic regime with respect to the Alfvén velocity taken in the background magnetic field.

Table 2

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), BT = δB + B0: (4)The electromagnetic fluctuating components δE and δB are provided by the MHD code. In this work, only protons are considered, hence m = mp and q = + e. Each particle is injected with a Lorentz factor γ0. This defines the particle’s Larmor radius rL0 = γ0mc2/eB0. We also define the particle’s synchrotron pulsation Ω0 = rL0/c. Equation (4) is normalized with respect to this initial value as (see Beresnyak et al. 2011, hereafter BYL11): We assume that the large-scale magnetic field B0 is aligned along ez. In the previous equations we have the following notations: , , . In practice, at least in the interstellar medium, the effect of the electric field over high-energy (multi TeV) CRs can be neglected and remains equal to 1.

When particles are submitted to radiative losses (or re-acceleration) 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 tpp ≃ (σppnHc)-1 for a medium of hydrogen density nH. We note that hereafter, as the particle energy is fixed, we denote rL0 and Ω0 to rL and Ω.

2.2.2. Integration schemes

We tested several integration schemes for Eqs. (5) and (6): a leap-frog scheme (by default implemented in RAMSES to treat the propagation of test particles in a gravitational field), a Runge-Kutta method of 5th order, and a Bulirsch-Stoer 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 leap-frog second-order 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 106 particles depends on the particle’s rigidity and the Alfvénic Mach number. It varies from 106 to a few 107 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 volume-averaged interpolation of the field strength. The volume can include the next eight neighbor points, this is the first-order cloud-in-cell (CIC) interpolation implemented by default in RAMSES or the next 64 neighbor points using a third-order 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 fast-magnetosonic modes are absent, the transport is controlled by the shear-Alfvén and the pseudo-Alfvén modes (the incompressible limit of slow-magnetosonic 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 back-reaction of CRs over the turbulence is considered in this work.

Beresnyak et al. (2011) have conducted 3D simulations at a resolution of 7683 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 Ma = 0.1 and in a regime of strong turbulence (or weak background field) with Ma = 1. We note that BYL11 calculate the Alfvén velocity with respect to the background magnetic field B0, hence Ma has to be interpreted as Ma0 = δB/B0 in our work. The injection scale of the turbulence in their simulation is (in box length units) Linj ≃ 0.2L. We note that in the sub-Alfvénic case, BYL11 use an elongated box with a main axis parallel to B0. Then, in that case Linj corresponds to the perpendicular turbulent outer-scale and the parallel outer-scale of the turbulence in the parallel direction is 10 Linj.

Beresnyak et al. (2011) have computed the variation of spatial diffusion coefficients with respect to the direction of the background magnetic field B0. 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 trans-Alfvénic turbulence, BYL11 find a parallel diffusion coefficient D independent of the particle’s rigidity at low rigidity ρ = rL/L ≤ 0.022 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 sub-Alfvénic and trans-Alfvé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 trans-Alfvénic case better. We find the same trend in the sub-Alfvé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 trans-Alfvénic cases.

Xu & Yan (2013) have investigated the propagation of CRs in compressible MHD turbulence and so the simulations retain the effects of shear-Alfvén modes and both magnetosonic modes. Here again the magnetostatic and test-particle limits have been retained.

Xu & Yan (2013) conducted 3D simulations at a resolution of 5123. 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 Ma = 0.19 to Ma = 0.73. XY13 have calculated the Alfvén velocity with respect to the total magnetic field BT hence their definition of Ma is similar to our. The injection scale of the turbulence in their simulation is (in box length units) Linj = 0.4 L.

Xu & Yan (2013) computed spatial diffusion coefficients with respect to the direction of the background magnetic field B0. 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 Ma. 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 >Linj 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 λ>Linj (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 Ma and ρ that extend the two above works significantly. We find power-law variations for λ and λ. The power-law 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 B0, x(t),y(t) determines the position of the particle in the directions perpendicular to the background magnetic field B0.

To proceed to mfp calculations, we have selected a set of Nr = 3 magnetic field realizations separated by at least two large-scale cascade times in order to ensure a statistical independence between two successive realizations. For each realization an ensemble of Np 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 ρ = rL/L = 0.037). We note that at low Ma the ballistic regime is present up to t ~ 103Ω-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 Np = 106 particles is used.

thumbnail Fig. 1

Time-dependent evolution of λ(t) (continuous lines) and λ(t) (dashed lines) for a normalized particle Larmor radius ρ = rL/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 Ma. 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 high-energy 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 EPeV and EGeV are the particle energies expressed in PeV and GeV units.

thumbnail 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 ρ = rL/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 Ma at all particle rigidities. The SoF solutions show faster variations with Ma than the CoF solutions. We find typical ratios λ(χ = 1) /λ(χ = 0) varying between ~100 at low Ma to ~1 as Ma → 1. We find that λ(χ = 1) /λ(χ = 0) varies from 4–5 for Ma< 0.6 to ~1 as Ma → 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 k0.097 deep in the inertial range of the turbulence. It corresponds to a scale reasonably far from the injection zone (typically k0.097L ~ 10,k0.097Linj ~ 2), as in our case Linj ~ 0.2L. Particles with rigidities in the inertial range are interesting for investigating the effect of resonant pitch-angle scattering over particle parallel mfps. The difficulty is that resonant interactions occurring at krL ~ 1 are important in the calculation of λ, and non-resonant interaction produced by the transit-time damping (TTD) mechanism due to fluctuations at scales larger than rL 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 Nr = 3 in all our calculations.

thumbnail Fig. 3

Upper figure: parallel mfp versus Alfvénic Mach numbers Ma at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: parallel mfp versus Alfvénic Mach numbers Ma 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 non-resonant pitch-angle 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)

thumbnail Fig. 4

Upper figure: perpendicular mfp versus Alfvénic Mach numbers Ma at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: perpendicular mfp versus Alfvénic Mach numbers Ma 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 Ma< 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 Linj. 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 Ma. 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 5123 (as can be clearly seen in Fig. 2). Wave-particle 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 Ma) have a uniform dependence on Ma. 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 δBB0. 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 slab-type turbulence with a steep inertial spectrum index5 a result obtained using a second order quasi-linear 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. Sub-Alfvénic turbulence is characterized by two regimes (Lazarian 2006): between Linj 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 kLinj ~ 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 λ/Linj> 1, then perturbations uncorrelated at scales larger than tr produce a perpendicular mfp . If conversely λ/Linj< 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 λ/Linj> 1 for curvilinear distances s along the magnetic field lines larger than a few tens of Linj. However, they described the existence of a super-diffusive perpendicular transport regime produced by the strong turbulence regime at intermediary scales s/Linj ~ 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 Ma 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 Ma points. This means that our results are consistent with QLT predictions in the low Ma regime and extend it in the high Ma 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 fast-magnetosonic) 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 Ma 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 Landau-synchrotron resonance (also known as gyro-resonance), 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 pitch-angle scattering. As reported above, we find typical ratios λ(χ = 1) /λ(χ = 0) ~ 10−100 at low Ma 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 gyro-resonance 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 fast-magnetosonic 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 gyro-resonance by Alfvén waves can contribute to the pitch-angle scattering substantially. Second, if the SoF-to-CoF mfp ratio is reduced by a factor of ~10, then the gyro-resonance 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 fast-magnetosonic 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 Ma< 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, λ/Linj varies in the range 1–103 (see the upper part of Fig. 3). At ρ = 0.01 we have λ/Linj> 100 and . So the result is not compatible with a scaling, but is still close to it. At such low rigidities and because λ/Linj ≫ 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 wave-particle 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 Ma. We do not find any clear trend of super-diffusive CR transport which could produce a larger final perpendicular mfps with respect to the diffusive regime. Super-diffusion is produced by the strong-Alfvénic turbulence (Lazarian & Yan 2014) in an interval which widens as Ma increases whereas perturbed scales associated with weak-Alfvénic turbulence which produces magnetic field line diffusion are restricted to scales ~Linj only. We think that we need a higher grid resolution at low Ma in order to include scales where GS spectrum develops in order to capture a possible effect of super-diffusion. 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 λ/Linj ~ 1, and not >1, so wave-particle 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 super-diffusion in that rigidity regime. At high Ma, 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)Linj. 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 fast-magnetosonic 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 Ma = 0.89 is a bit under (above) the extrapolation of in the CoF case (SoF case). At high Ma indeed the use of QLT is questionable.

thumbnail Fig. 5

Upper figure: parallel mean free path dependence with Ma at ρ = 0.097 for χ = 0,0.5,1. Lower figure: perpendicular mean free path dependence with Ma at ρ = 0.097 for χ = 0,0.5,1.

Open with DEXTER

thumbnail 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 Ma< 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 (fast-magnetosonic) waves preferentially, our results can give some qualitative hints about the respective effect of gyro-resonance 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 gyro-resonance by the Alfvénic turbulence. Also, the ratio of the parallel SoF and CoF mfps decreases, which can be due to the effect of gyro-resonance of Alfvén and fast-magnetosonic 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 λ/Linj< 1 or λ/Linj> 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 super-diffusion 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 Ma ~ 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 Ma at ρ = 0.097 for χ = 0,0.5, and 1.

It can be seen that in both cases (parallel and perpendicular mfps), especially for Ma ≤ 0.5, the MoF solutions at χ = 0.5 are closer to the SoF solutions. This is supported by the value of rc and the fraction of compressible modes (see Sect. 2.1.1): in the MoF case rc = 1/3 is closer to rc = 0 obtained for χ = 1 than rc = 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 ρ = rL/L.

thumbnail Fig. 7

Upper figure: parallel mfps versus particle normalized Larmor radii ρ at different Ma in the SoF case (χ = 1). Insets: slopes obtained for different indices. Lower figure: parallel mfp versus particle normalized Larmor radius ρ at different Ma in the CoF case (χ = 0).

Open with DEXTER

Table 3

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 ρ = rL/L.

thumbnail Fig. 8

Upper figure: perpendicular mfp versus particle normalized Larmor radius ρ at different Ma values in the SoF case (χ = 1). Lower figure: perpendicular mfp versus particle normalized Larmor radius ρ at different Ma values in the CoF (χ = 0).

Open with DEXTER

Table 4

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 Ma 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 Ma (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 = Pg/Pm ~ (Ma/Ms)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 Ma = 1 and Ma = 0.1, but presented their results for λ at Ma = 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 (7683 instead of 5123). Again, we think these results should be considered with some caution especially for λ. For λ, λρ1/3 at Ma = 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 Ma = 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 Goldreich-Sridhar (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/B0 = 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 Ma< 18. Yan & Lazarian (2008) instead reported on a dominance of fast-magnetosonic 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 non-linear 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 = B0 (even if the fit is restricted to almost one decade in rigidity). At lower Ma, only if λ<LinjYan & 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 Ma = 0.89 which is marginally compatible. The mfp index in the regime ρ< 0.1 shows some hardening from low to high Ma, 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 Ma, 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 Ma ≤ 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 (Alouani-Bibi & le Roux 2014). However, this possibility is questionable since we did not detect any subdiffusive paths at low rigidities and low Ma. Inverted indices can be associated with a propagation of CRs in a damped turbulence (Yan & Lazarian 2008; Shalchi & Büsching 2010). Above Ma = 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 fast-magnetosonic 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 Ma = 0.34 and Ma = 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 Ma.

thumbnail Fig. 9

Magnetic energy density power spectra for Ma = 0.88 (blue continuous line) and Ma = 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 Ma = 0.34 and Ma = 0.89. We checked that the stationarity was reached. If the spectrum at Ma = 0.89 has an index s ~ 1.5 the spectrum at Ma = 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 Ma = 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 (Ma,Ms) 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 sub-Alfvénic turbulence (see Sect. 3.3). Hence, if Ma = 0.3, weak sub-Alfvé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.

thumbnail Fig. 10

Magnetic energy density power spectra for Ma = 0.88 (blue continuous line) and Ma = 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 high-rigidity 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 Linj ~ 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 rL ~ 0.5 Linj as Linj = 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 Ma). At low Ma, the index is also compatible with the solution obtained at level 9 (see Table 4). At Ma = 0.66 and Ma = 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 Ma ~ 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 Ma, 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 counter-propagating 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 Ma) but also with non-linear models of CR diffusion in fast-magnetosonic 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 Ma 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 sub-Alfvé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 Ma = 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 high-rigidity points where the contribution of pitch-angle 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 Ma = 0.89. If χ = 1 and Ma< 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 sub-Alfvé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.

thumbnail 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 + (λ/rL)2).

Open with DEXTER

4. Summary and conclusions

In this work we have developed a series of kinetic-MHD 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 particle-in-cell 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 Ma dependences with respect to the results of XY13, but our results are compatible at Ma> 0.7. Even if we found a similar trend to XY13 for the dependence of λ on Ma, 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 λ/Linj is neither ≫ 1 nor ≪ 1 in our solutions. We did not find a clear signature of super-diffusion in the CR perpendicular mfp produced by GS-type 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 fast-magnetosonic 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 sub-Alfvé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 Ma. 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 small-scale 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 Ma. 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 non-linear 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 gyro-resonance. 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 Ma ≤ 0.1 Alfvénic Mach number regime. It would be also interesting to perform high-resolution 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 super-Alfvénic (Ma ≥ 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 large-scale turbulence which is, for instance the case of high-energy CRs (TeV and beyond) in the ISM. The transport of low-energy particles in the ISM requires adding a component of self-generated waves (see the discussion in Sect. 6 of XY13). As the pressure imparted in the low-energy CR is as important as the gas and magnetic pressure, this would be possible using PIC-MHD 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.


1

and k are the perturbed wave numbers defined along the local magnetic field BL and the total magnetic field BT, respectively. We also define the global or background magnetic B0, 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.

2

We adopt the following notations: is the reduced rigidity with respect to the turbulence outer scale and ρ is the reduced rigidity in cube units; in BYL11 and in our work we have .

3

Jokipii & Parker (1968) instead predicted and .

4

Casse et al. (2002) defined η as a quantity which corresponds to η2 in our text.

5

This is the 1D wave number index, denoted s in Sect. 3.5.4. Deviations from QLT are obtained for s> 2.

6

Yan & Lazarian (2008) used Ma0 instead of Ma.

7

βp is the plasma parameter, Pp and Pm are the gas thermal pressure and the magnetic pressure respectively.

8

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 super-diffusive 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

All Tables

Table 1

MHD simulations used in this work.

Table 2

Correspondence between the parameters η and η0.

Table 3

Parallel mean free path indices and errors for SoF and CoF geometries.

Table 4

Perpendicular mean free path indices and errors for both forcing geometries.

All Figures

thumbnail Fig. 1

Time-dependent evolution of λ(t) (continuous lines) and λ(t) (dashed lines) for a normalized particle Larmor radius ρ = rL/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
thumbnail 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 ρ = rL/L.

Open with DEXTER
In the text
thumbnail Fig. 3

Upper figure: parallel mfp versus Alfvénic Mach numbers Ma at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: parallel mfp versus Alfvénic Mach numbers Ma 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
thumbnail Fig. 4

Upper figure: perpendicular mfp versus Alfvénic Mach numbers Ma at different particle rigidities in the case of a solenoidal forcing (χ = 1). Lower figure: perpendicular mfp versus Alfvénic Mach numbers Ma at different particle rigidities in the case of a compressible forcing (χ = 0).

Open with DEXTER
In the text
thumbnail Fig. 5

Upper figure: parallel mean free path dependence with Ma at ρ = 0.097 for χ = 0,0.5,1. Lower figure: perpendicular mean free path dependence with Ma at ρ = 0.097 for χ = 0,0.5,1.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 7

Upper figure: parallel mfps versus particle normalized Larmor radii ρ at different Ma in the SoF case (χ = 1). Insets: slopes obtained for different indices. Lower figure: parallel mfp versus particle normalized Larmor radius ρ at different Ma in the CoF case (χ = 0).

Open with DEXTER
In the text
thumbnail Fig. 8

Upper figure: perpendicular mfp versus particle normalized Larmor radius ρ at different Ma values in the SoF case (χ = 1). Lower figure: perpendicular mfp versus particle normalized Larmor radius ρ at different Ma values in the CoF (χ = 0).

Open with DEXTER
In the text
thumbnail Fig. 9

Magnetic energy density power spectra for Ma = 0.88 (blue continuous line) and Ma = 0.34 (red continuous line) in the CoF case.

Open with DEXTER
In the text
thumbnail Fig. 10

Magnetic energy density power spectra for Ma = 0.88 (blue continuous line) and Ma = 0.34 (red continuous line) in the SoF case.

Open with DEXTER
In the text
thumbnail 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 + (λ/rL)2).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.