A NUMERICAL ASSESSMENT OF COSMIC-RAY ENERGY DIFFUSION THROUGH TURBULENT MEDIA

Until recently, the widely held paradigm for the origin of cosmic rays, at least below the knee at roughly 1 PeV, promoted the view that all of the intragalactic injection occurred via first-order Fermi acceleration in supernova shells. However, any direct evidence for this view is meager and equivocal. More recently, data from balloon-borne experiments have refuted the expectation from supernova acceleration schemes that the cosmic-ray spectrum ought to be structureless and universal (Wefel 1988). The latest measurements with PAMELA confirm and extend the balloon-based claims (Adriani et al. 2011). These data seem to call for a diverse variety of acceleration sites and mechanisms throughout the Galaxy.

In recent work (Melia & Fatuzzo 2011; Fatuzzo & Melia 2012a, 2012b), we have begun to assess the feasibility of stochastic acceleration within turbulent magnetized regions using highly detailed simulations of individual particle trajectories. A principal goal of this work is to accurately determine the spatial and energy diffusion coefficients of cosmic-ray protons in a broad range of environments, e.g., inside molecular clouds and the more tenuous intercloud medium. The spatial and energy diffusion coefficients calculated over a broad range of parameter space may be used, e.g., to compare estimates of the time required to energize protons up to TeV energies with the escape and cooling times throughout the interstellar medium (ISM). In previous applications, this approach has allowed us to conclude that protons in the intercloud medium at the galactic center can be energized up to the 1–10 TeV energies required to account for the observed HESS emission in this region (Aharonian et al. 2006; Liu et al. 2006a; Ballantyne et al. 2007; Wommer et al. 2008; see also Markoff et al. 1997, 1999; Liu & Melia 2001). Stochastic particle acceleration by magnetic turbulence appears to be a viable mechanism for cosmic-ray production at the galactic center (Liu et al. 2006b).

Developing an understanding of how the spatial diffusion coefficients depend on the physical environment and energy of the particles was the subject of our previous work (Fatuzzo et al. 2010). The focus of the present paper is specifically to determine the diffusion of cosmic rays in energy space, with its attendant dependence on all of the critical physical characteristics of the medium, such as the magnetic intensity, degree of turbulence, and size of the fluctuations. As before, we do this by using a modified numerically based formalism developed for the general study of cosmic-ray diffusion by Giacalone & Jokipii (1994). This approach has already been used successfully in several other contexts (see, e.g., O'Sullivan et al. 2009; Fatuzzo & Melia 2012a, 2012b). Here, we will extend this robust numerically based framework for the general analysis of stochastic acceleration of cosmic-rays by the turbulent electric fields generated along with the time-dependent turbulent magnetic field in a dynamically active medium.

We stress that our adopted model of turbulence as a superposition of Alfvénic like waves with linear dispersion relations is not meant to fully describe MHD turbulence in the ISM. For example, the model does not account for the fact that magnetic fluctuations decorrelate due to non-linear interactions before they can propagate over distances of multiple wavelengths (Goldreich & Sridhar 1995)—an effect that leads to resonance broadening—and as such, influences how thermal particles interact with turbulence (Lynn et al. 2012, 2013). However, our intent is to adopt a formalism that adequately describes turbulence as seen locally by highly relativistic cosmic rays. Since relativistic particles have speeds that are much greater than the Alfvén speeds considered in our analysis (which is then the limit in which our results are expected to be valid), they should not be sensitive to dynamical processes that occur on MHD timescales. In essence, we are relying on basic principles (e.g., the scaling laws between intensity and wavelength) to capture the global features of MHD turbulence that affect the propagation and acceleration of high energy cosmic rays.

The work presented here extends the results of previous works, most notably, that of O'Sullivan et al. (2009), as it significantly broadens the explored parameter space and also considers anisotropically distributed wave vectors (Goldreich & Sridhar 1995; Cho & Vishniac 2000). We focus primarily on strong turbulence (δB ∼ B) for which quasilinear theory is not applicable, although we do present results for isotropic weak turbulence as a consistency check. Our paper is organized as follows. The pertinent properties of the medium through which the particles diffuse are outlined in Section 2. The scheme for generating the turbulent magnetic and electric fields is presented in Section 3, along with the equations that govern the motion of cosmic rays. The basic elements of stochastic acceleration are discussed in Section 4, and the results of our work are presented in Section 5. The conclusions of our work are summarized in Section 6.

The physical parameters found throughout the intergalactic and interstellar media have values that span several orders of magnitude. Of particular interest to our study, the most vacuous regions of the intergalactic medium have particle number densities n ≲ 10⁻³ cm⁻³ and magnetic field strengths B ≲ 0.1 μG (see, e.g., Kronberg 1994; Fraschetti & Melia 2008). In contrast, the denser regions near the supermassive black hole at the galactic center have densities n ≳ 10¹² cm⁻³ and field strengths B ≳ 1 G (Ruffert & Melia 1994; Falcke & Melia 1997; Kowalenko & Melia 1999; Melia 2007; see also Misra & Melia 1993 for the case of stellar-mass size black holes).

Exactly how the magnetic field is partitioned within these various media is not yet known, but there is a general relation between density and field strength. In the simplest case where flux freezing applies (say in the ISM), the magnetic field strength B would scale with gas density n according to B∝n^1/2. It is noteworthy, then, that an analysis of magnetic field strengths measured in molecular clouds yields a relation between B and n of the form

$$\begin{equation} B \sim 10 \,\mu {\rm G} \left({n \over 10^2\, {\rm cm}^{-3}} \right)^{0.47}, \end{equation} \tag{ 1 }$$

though with a significant amount of scatter in the data used to produce this fit (Crutcher 1999; see also Fatuzzo et al. 2006 and references therein). This result is consistent with the idea that nonthermal linewidths, measured to be ∼1 km s⁻¹ throughout the cloud environment (e.g., Lada et al. 1991), arise from MHD fluctuations.

Of course, even among molecular clouds, the physical environment can be quite different depending on location. Near the galactic center, the molecular clouds are considerably different from their counterparts in the disk. For example, the average molecular hydrogen number density over the Sgr B complex—the largest molecular cloud complex (Lis & Goldsmith 1989; Lis & Goldsmith 1990; Paglione et al. 1998) near the galactic center—has a density 3–10 × 10³ cm⁻³. Like its traditional counterparts, Sgr B displays a highly nonlinear structure, containing two bright sub-regions, Sgr B1 and Sgr B2, the latter having an average molecular density of ∼10⁶ cm⁻³, and containing three dense (n ∼ 10^7.3–8 cm⁻³), small (r ∼ 0.1 pc) cores—labeled north, main, and south. These cores also show considerable structure, containing numerous ultra-compact and hyper-compact H ii regions. As such, the densities associated with the galactic-center molecular clouds are about two orders of magnitude greater than those in the disk of our galaxy.

The exact nature of magnetic turbulence in these environments is itself not well constrained, though magnetic fluctuations typically have a power-law spectrum. Their intensity at a given wavelength scales according to (δB_λ)² ∼ λ^{Γ − 1}, with values of Γ = 1 (Bohm), Γ = 3/2 (Kraichnan), or Γ = 5/3 (Kolmogorov) often adopted. In addition, the range in wavelengths over which these fluctuations occurs is not well known, though it is reasonable to assume that the upper end corresponds to the lengthscale over which the fluctuations are generated. (For example, in the ISM, the turbulence is generated by supernova remnants and stellar-wind collisions, so one might expect the longest wavelength to be on the order of several parsecs or less; see, e.g., Coker & Melia 1997; Melia & Coker 1999.) The lower end must be smaller than the characteristic length scale associated with the particle motion (i.e., the gyration radius), under the assumption that magnetic energy ultimately dissipates into plasma energy via its coupling to these particles.

We explore how cosmic rays diffuse through a homogeneous hydrogen gas of mass density ρ = nm_p threaded by a uniform static background field B₀, on which magnetic and electric fluctuations propagate. From our numerical simulations, we then compute energy diffusion coefficients covering a broad range of particle energies and parameter space expected to span the great diversity of environments observed both within our galaxy and in the intergalactic medium. We focus primarily on strong turbulence, for which the energy density of the turbulent fields is comparable to that of the background fields.

A standard numerical approach to studying cosmic-ray diffusion treats the spatially fluctuating magnetic component δB as the superposition of a large number of randomly polarized transverse waves with wavelengths λ_n = 2π/k_n, logarithmically spaced between λ_min and λ_max (e.g., Giacalone & Jokipii 1994; Casse et al. 2002; Fatuzzo et al. 2010). Adopting a static turbulent field removes the necessity of specifying a dispersion relation between the wavevectors k_n and their corresponding angular frequencies ω_n. This approach appears suitable for considering highly non-linear turbulence (δB ≫ B₀), or simply an environment without a background component. Of course, turbulent magnetic fields in cosmic environments are not static. Nevertheless, a static formalism in spatial diffusion calculations of relativistic particles seems justified for environments in which the Alfvén speed is much smaller than the speed of light.

The situation in this paper is quite different: we are focusing on the energy diffusion of cosmic rays propagating through a turbulent magnetic field, which requires the use of a time-dependent formalism in order to self-consistently include the fluctuating electric fields that must also be present (say, from Faraday's law). At present, such a theory of MHD turbulence in the ISM remains elusive. Nevertheless, it is generally understood that turbulence is driven from a cascade of longer wavelengths to shorter wavelengths as a result of wave–wave interactions. For strong MHD turbulence in a uniform medium, this cascade seemingly produces eddies on small spatial scales that are elongated in the direction of the underlying magnetic field, so that the components of the wave vector along (k_||) and across (k_⊥) the underlying field direction are related by the expression $k_{||}\propto k_\perp ^{2/3}$, with a Kolmogorov energy spectrum that scales as $k_\perp ^{-5/3}$ (Goldreich & Sridhar 1995; Cho & Vishniac 2000).

It is beyond the scope of this paper to develop a self-consistent theory of MHD turbulence in the ISM that includes electric field fluctuation. We therefore adopt the formalism of O'Sullivan et al. (2009). Specifically, we assume a medium represented by a nonviscous, perfectly conducting fluid threaded by a uniform static field B₀, and use linear MHD theory as a guide. In the linear regime, one can encounter three types of MHD waves: Alfvén, fast and slow. As in O'Sullivan et al. (2009), we here consider only Alfvén waves, for which the turbulent magnetic field may be written as a sum of N randomly directed waves

$$\begin{equation} \delta {\bf B} = \sum _{n=1}^N {\bf A_n} \, e^{i({\bf k_n} \cdot {\bf r}- \omega _n t+\beta _n)}\,. \end{equation} \tag{ 2 }$$

As noted already, this formalism does not adequately describe decorrelation effects known to be important for thermal cosmic rays (Lynn et al. 2012, 2013). However, it should serve as a reasonable model for the turbulence experienced by particles with velocities much greater than the Alfvén speed v_A, as such particles would be expected to travel over many correlation lengths (∼0.1λ_max) in an Alfvén time τ_A ∼ λ_max/v_A.

To keep the analysis as broad as possible, we focus on an isotropic turbulence spectrum, but we also perform a suite of experiments with anisotropic turbulence as informed by the results of Goldreich & Sridhar (1995). For the isotropic case, the direction of each propagation vector k_n is set through a random choice of polar angles θ_n and ϕ_n, and the phase of each term is set through a random choice of β_n. The sum has N = N_klog₁₀[k_max/k_min] terms, with k_n evenly spaced on a logarithmic scale between k_min = 2π/λ_max and k_max = 2π/λ_min. A value of N_k = 25 appears to provide enough terms to suitably model what is really a continuous rather than a discrete system. The appropriate choice of Γ in the scaling

$$\begin{equation} A_n^2 = A_1 ^2\left[{k_n \over k_1} \right]^{-\Gamma } {\Delta k_n\over \Delta k_1} = A_1^2\left[{k_n \over k_1} \right]^{-\Gamma +1} \end{equation} \tag{ 3 }$$

sets the desired spectrum of magnetic turbulence (e.g., Γ = 3/2 for Kraichnan and 5/3 for Kolmogorov turbulence). For our adopted scheme, the value of Δk_n/k_n is the same for all n. The normalization for A₁ is set by the parameter

$$\begin{equation} \eta = {\langle \delta B^2 \rangle \over B_0^2}, \end{equation} \tag{ 4 }$$

defined as the ratio of magnetic energy density in the turbulent component to that of the static background field B₀.

For the anisotropic case, the direction of each perpendicular wavevector k_n⊥ is set through a random choice of azimuthal angle ϕ_n, and the phase of each term is again set through a random choice of β_n. The corresponding parallel component of the wavevector is then set through the relation

$$\begin{equation} k_{n||} = \pm {\sqrt{2}\over 2} k_{1\perp }^{1/3}\,k_{n\perp }^{2/3}\, \end{equation} \tag{ 5 }$$

(with a randomly chosen sign), such that

$$\begin{equation} k_1 = \sqrt{k_{1\perp }^2+k_{1||}^2} = {2\pi \over \lambda _{{\rm max}}}. \end{equation} \tag{ 6 }$$

Consistent with the results of Goldreich & Sridhar (1995), we only consider a Kolmogorov profile, so that

$$\begin{equation} A_n^2 = A_1^2\left[{k_{n\perp } \over k_{1\perp }} \right]^{-2/3}\,. \end{equation} \tag{ 7 }$$

Since Alfvén waves don't compress the fluid through which they propagate, their fluid velocity v satisfies the condition k · v = 0. In addition, v · B₀ = 0 for Alfvén waves. The fluid velocity associated with the nth term in Equation (2) is therefore

$$\begin{equation} \delta {\bf v_n} = \pm A_n \, {v_A\over B_0} \, {{\bf B_0} \times {\bf k_n} \over | {\bf B_0} \times {\bf k_n} |} \, e^{i({\bf k_n} \cdot {\bf r}-\omega _n t+\beta _n)}\,, \end{equation} \tag{ 8 }$$

where the sign is chosen randomly for each term in the sum. The dispersion relation for Alfveńic waves is given by the expression

$$\begin{equation} \omega _n = v_A |k_{n\Vert }|\,. \end{equation} \tag{ 9 }$$

Each wave has a magnetic field given by the linear form of Ampère's law,

$$\begin{equation} {\bf A_n} = \mp A_n \, {{\bf k_n} \cdot {\bf B_0} \over | {\bf k_n} \cdot {\bf B_0} |} \, {{\bf B_0} \times {\bf k_n} \over | {\bf B_0} \times {\bf k_n} |} \,, \end{equation} \tag{ 10 }$$

which is identical to that of O'Sullivan et al. (2009).

Insofar as the electric fields are concerned, if we were to naively extrapolate from the results of linear MHD theory, the total electric field δE associated with the turbulent magnetic field in Equation (2) would be given by a sum over the terms

$$\begin{equation} \delta {\bf E_n}= -\delta {\bf v_n} \times {\bf B_0}\,. \end{equation} \tag{ 11 }$$

Notice that δE · B₀ = 0, but the second order term δE · δB ≠ 0. The presence of an electric field component parallel to the magnetic field in this second order term can significantly increase the acceleration efficiency artificially, especially if the formalism is extended to the nonlinear regime (δB ∼B₀). However, the ISM is highly conductive, so any electric field component parallel to the magnetic field should be quickly quenched. O'Sullivan et al. (2009) circumvented this problem by first obtaining the total fluid velocity δv via the summation

$$\begin{equation} \delta {\bf v} = \sum _{n=1}^N \delta {\bf v_n}\,, \end{equation} \tag{ 12 }$$

and then using the MHD condition to set the total electric field:

$$\begin{equation} \delta {\bf E} = - {\delta {\bf v}\over c} \times {\bf B}, \end{equation} \tag{ 13 }$$

where B = B₀ + δB. This is the procedure we too will use here.

With these electric and magnetic field components, one may then solve the Lorentz force equation

$$\begin{equation} {d\over dt} (\gamma m_p {\bf v})= e \left[\delta {\bf E} + {{\bf v} \times {\bf B}\over c}\right], \end{equation} \tag{ 14 }$$

with

$$\begin{equation} {d {\bf r}\over dt} = \bf v, \end{equation} \tag{ 15 }$$

to determine the motion of a relativistic charged proton with Lorentz factor γ through the turbulent medium. The solutions to these equations are not sensitive to the value of λ_min so long as the particle's gyration radius R_g = γmc²/(eB) ≫ λ_min (Fatuzzo et al. 2010). As such, we set λ_min = 0.1γmc²/(eB₀) in all our simulations. To completely specify the physical state of an environment to be studied, we must therefore provide values for the parameters B₀, n, λ_max, Γ and η.

The motion of a charged particle in perpendicular uniform electric and magnetic fields represents a fundamental topic in electrodynamics and is well understood. It is therefore easy to show that in general, the charge gains and loses kinetic energy in cyclic fashion as it "drifts" in the E × B direction.

In the turbulent medium in which the magnetic and electric fields are perpendicular, a single particle's energy, characterized by Δγ/γ₀ = γ/γ₀ − 1, therefore exhibits a random-walk like behavior, as shown in Figure 1. The energy distribution of an ensemble of particles injected with the same Lorentz factor γ₀ thus broadens as the particles sample the turbulent nature of the accelerating electric fields. Interestingly, these distributions appear Gaussian for all of the isotropic and η = 1 anisotropic cases explored in our analysis. This point is illustrated by Figure 2, which shows the distribution of Δγ values at time t = 1000λ_max/c for an ensemble of 1000 particles injected with γ₀ = 10⁵ into a medium defined by the same parameters used to calculate the energy evolution shown in Figure 1, but with an anisotropic turbulence profile. In such cases, one can therefore quantify the stochastic acceleration of particles in turbulent fields through the dispersion of the resulting distributions of initially mono-energetic particles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast, the particle distributions obtained for anisotropic, weak turbulence (η ≪ 1) cases have wings that are much broader than a Gaussian distribution, so that their dispersion significantly overestimates the true distribution width. This point is illustrated in Figure 3, which shows the distribution of Δγ values at time t = 1000λ_max/c for an ensemble of 1000 particles injected with γ₀ = 10⁵ into a medium defined by the same parameters used to calculate the distribution shown in Figure 2, but with η = 0.01.

Zoom In Zoom Out Reset image size

To further compare strong and weak isotropic and anisotropic turbulence, we plot in Figure 4 the dispersion $\sigma _\gamma \equiv \sqrt{\langle \Delta \gamma ^2\rangle }$ of the Δγ distribution as a function of time for the particles used to generate Figures 2 and 3 (for which the turbulence was anisotropic) along with their isotropic counterparts. As found throughout our analysis, there is little difference between the isotropic and anisotropic cases when the turbulence is strong (η = 1). However, anisotropic turbulence appears to be significantly less effective at energizing particles than isotropic turbulence when η ≪ 1, in agreement with results obtained in the quasi-linear approximation (Chandran 2000; see also Yan & Lazarian 2002).

Zoom In Zoom Out Reset image size

Clearly, the chaotic nature of motion through turbulent fields necessitates a statistical analysis. We define a single experiment as a numerical investigation of particle dynamics through a given environment (as defined by the parameters B₀, n, Γ, λ_max, and η) over a broad range of particle injection energies (as defined by the Lorentz factor γ₀). As expected from the random nature of stochastic acceleration, $\sigma _\gamma \propto \sqrt{t}$ once particles have had a chance to sample the turbulent nature of the underlying electric fields, i.e., for t ≳ λ_max/c. This in turn means that the energy diffusion coefficient D_γ ≡ 〈Δγ²〉/(2Δt) can be calculated by using an integration time Δt ≳ λ_max/c. For each particle energy, we numerically integrate the equations of motion for a time Δt = 10λ_max/c for N_p = 1000 protons randomly injected from the origin. Each particle samples its own unique magnetic field structure (i.e., the values of β_n, θ_n, ϕ_n, and the choice of a ± are chosen randomly for each particle for the isotropic case).

We use the procedure described above to carry out a suite of experiments sampling a broad region of parameter space for isotropic turbulence, and perform a limited complimentary set of experiments for anisotropic, strong (η = 1) Kolmogorov turbulence. We limit our analysis to particle energies for which the gyration radius R_g falls comfortably below the maximum turbulence wavelength λ_max, so that particles actually diffuse through the turbulent medium. A principal goal of this paper is to determine the relationship between the energy diffusion coefficient D_γ and the particle's energy for each experiment. We therefore fit the numerical data using a power law

$$\begin{equation} D_\gamma = D_{\gamma 0} \gamma ^{\alpha }\,. \end{equation} \tag{ 16 }$$

The parameters for each experiment and corresponding values of D_γ0 and α obtained through the fits are summarized in Table 1. As a reference, we then compare these empirical fits to the quasi-linear expression

$$\begin{equation} D^{ql}_\gamma \approx {v_A^2\over c^2} \left({\delta B\over B_0}\right)^2 \left({R_g\over \lambda _{{\rm max}}}\right)^{\Gamma -1} {\gamma ^2 c\over R_g} \,, \end{equation} \tag{ 17 }$$

(Schlickeiser 1989, O'Sullivan et al. 2009). We note that in terms of our parameters, this predicted expression reduces to the simpler form

$$\begin{equation} D^{ql}_\gamma \propto n^{-1} B_0^{4-\Gamma } \eta \,\lambda _{{\rm max}}^{1-\Gamma }\, \gamma ^\Gamma \,. \end{equation} \tag{ 18 }$$

We plot the results of Experiments 1–4 in Figure 5 and Experiments 1, 4, and 5–8 in Figures 6 and 7 to illustrate how the diffusion coefficient index changes between the strong (η ∼ 1) and weak (η ≪ 1) turbulence limits. As predicted by quasi-linear theory, α ≈ Γ when η ≪ 1 for isotropic turbulence, with the best match occurring for Kolmogorov diffusion. However, Equation (17) overestimates the value of D_γ0 by as much as two orders of magnitude. In addition, the indices decrease as the turbulence grows stronger, and quasi-linear theory does not appear to adequately describe the energy dependence for Kolmogorov (Γ = 5/3) and Kraichnan (Γ = 3/2) diffusion when η ∼ 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next consider how the diffusion coefficients for Kolmogorov (Γ = 5/3), Kraichnan (Γ = 3/2) and Bohm (Γ = 1) turbulence scale with n, B₀ and λ_max, but limit our analysis to the strong turbulence limit (η = 1). Our results for Kolmogorov turbulence (derived from Experiments 1 and 9–20) are shown in Figures 8–13. Specifically, Figure 8 illustrates the energy dependence of D_γ for four values of ambient density n. As expected, a greater density yields a smaller Alfvén speed, and therefore a reduced energy diffusion. However, the index α is not sensitive to n. We note that Equation (17) overestimates the energy diffusion coefficient for γ₀ ≳ 10³, and underestimates the energy diffusion coefficient for γ₀ ≲ 10³. Figure 9 then illustrates the relation between D_γ and n for the value γ₀ = 10⁵, where the fits from Figure 8 (along with fits to the data not shown in that figure) are used to obtain the plotted data points. The relationship between D_γ and n clearly takes on the same scaling obtained from quasi-linear theory, i.e., D_γ∝n⁻¹. This result is expected since D_γ∝(Δγ)²∝ (δE)²∝ $v_A^2 \propto n^{-1}$. We therefore assume that this scaling holds for all values of Γ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 illustrates the energy dependence of D_γ for five values of background field strengths B₀. Since a greater magnetic field strength leads to a greater electric field strength, energy diffusion increases with field strength. Again, the index α is not sensitive to B₀, but we find that Equation (17) overestimates the energy diffusion for larger values of γ₀, and underestimates the energy diffusion for smaller values of γ₀, although the transition varies depending on the field strength. Figure 11 illustrates the relation between D_γ and B₀ for the value γ₀ = 10⁵, where the fits from Figure 10 are used to obtain the plotted data points. We note that the scaling $D_\gamma \propto B_0^{2.50}$ obtained from our results differs from quasi-linear theory ($D_\gamma \propto B_0^{2.33}$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 illustrates the energy dependence of D_γ for three values of λ_max. Again, the index α is not sensitive to λ_max, and once again, Equation (17) overestimates the energy diffusion for larger values of γ₀, and underestimates the energy diffusion for smaller values of γ₀. Figure 13 illustrates the relation between D_γ and λ_max for the value γ₀ = 10⁵, where the fits from Figure 12 are used to obtain the plotted data points. As was the case with the field, the scaling obtained ($D_\gamma \propto \lambda _{{\rm max}}^{-0.47}$) differs from that predicted by quasi-linear theory ($D_\gamma \propto \lambda _{{\rm max}}^{-0.67}$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Corresponding results for isotropic Kraichnan turbulence (derived from Experiments 5 and 21–26) are shown in Figures 14–17, and for isotropic Bohm turbulence (derived from Experiments 7 and 27–32) are shown in Figures 18–21. The analysis for these cases mirrors that described above for Kolmogorov turbulence, and the same general results are obtained.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our results indicate that simple scaling laws of the form D_γ = D_γ0 γ^α, where

$$\begin{equation} D_{\gamma 0}= D_0 \left({n\over 1\,{\rm cm}^{-3}}\right)^{-1} \left({B_0\over 1\,\mu {\rm G}}\right)^\delta \left({\lambda _{{\rm max}}\over 1 \,{\rm pc}}\right)^\kappa \end{equation} \tag{ 19 }$$

can be used to obtain values of the energy diffusion coefficient over a wide range of parameters that pertain to turbulent interstellar and intergalactic environments, so long as the particle gyration radius λ_min ≪ R_g ≪ λ_max. In addition, the same expressions can be used to describe the isotropic and anisotropic cases in the limit of strong turbulence (η ≈ 1).

We obtain values of D₀, α, δ, and κ in the strong turbulent limit (η = 1) for the three turbulence profiles considered in our work by using the results presented above. Specifically, final values of D₀ and α, along with their 1σ errors, are obtained by finding the mean and standard deviations of the D₀ and α values from the corresponding experiments listed in Table 1. The values of δ and κ are obtained from the fits to the data shown in Figure 11 and Figure 13 for Γ = 5/3, Figure 15 and Figure 17 for Γ = 3/2, and Figure 19 and Figure 21 for Γ = 1. The results are summarized in Table 2.

We conclude our analysis by applying the results of this work toward obtaining estimates of the acceleration time τ_acc ≡ γ²/D_γ required to energize protons up to energies of 1 TeV in molecular cloud environments. To keep this analysis as simple as possible, we assume that the magnetic field scales with density as given by Equation (1). We also assume that the maximum turbulence wavelength λ_max scales as the cloud size. Following Larson's law (Larson 1981), we obtain the relation

$$\begin{equation} \lambda _{{\rm max}} = 1600 \,{\rm pc} \left({n_{H_2}\over 1 \,{\rm cm}^{-3}}\right)^{-0.91}\,. \end{equation} \tag{ 20 }$$

A plot of the acceleration time τ_acc as a function of molecular hydrogen density $n_{H_2}$ is shown in Figure 22 for Kolmogorov, Kraichnan, and Bohm turbulence. Given that clouds are not expected to last for more than ∼10 Myr, TeV cosmic ray production from stochastic acceleration of turbulent magnetic fields can clearly be ruled our for normal molecular cloud environments. We note, however, that the molecular clouds near the galactic center (GC) have fairly extreme environments. As shown in Fatuzzo & Melia (2012b), the acceleration times is the GC environment are considerably shorter, as illustrated by the solid circle (inter cloud region at the GC) and solid square (molecular cloud at the GC) shown in Figure 22.

Zoom In Zoom Out Reset image size

We have used a detailed numerical simulation to determine the spatial and spectral profiles of cosmic rays diffusing with arbitrary energy through turbulent media characterized by a broad range of magnetic fields, turbulence strength, fluctuation size, and ambient particle density, for both isotropic and anisotropic turbulence. This study has expanded considerably from the initial attempt at simulating the behavior of relativistic particle propagation through the galactic-center environment, where they encounter a variety of physical conditions, inside and outside of the molecular gas in that region. Our goal throughout this exercise has been to avoid using "standard" techniques, e.g., quasi-linear theory or the diffusion equation, all of which are often subject to unknown factors that delimit the applicability of these approaches to real systems (but see also Nayakshin & Melia 1998; Wolfe & Melia 2006). Indeed, one of the principal benefits of the technique we have developed in this work is an accurate determination of the spatial and energy diffusion coefficients that in turn may be used in these other approaches without the need to guess or estimate their normalization and energy dependence.

As the sensitivity and spectral range of high-energy observatories continue to improve, the need to accurately simulate the propagation of relativistic particles through turbulent media arises in an ever increasing range of environments, from the ISM, to compact accretion regions surrounding supermassive black holes (such as Sgr A* at the galactic center), to the hot intracluster gas, and in even more exotic structures, such as the Fermi bubbles straddling the center of the Milky Way (Crocker & Aharonian 2011). However, the physical conditions characterizing these regions change considerably from one environment to the next. For example, the magnetic intensity may be as large as several Gauss near accreting black holes, but smaller than 0.01 μG in the intergalactic medium. We now know that quasi-linear theory is only approximately valid even in weakly turbulent environments, let alone regions where the turbulent magnetic energy is comparable to, or bigger than, the underlying uniform component. Worse, it is often difficult to estimate the absolute value of the diffusion coefficients without resorting to observational data which, however, are sometimes difficult to get (as in the case of the giant radio lobes in FR II galaxies).

These are among the reasons we have embarked on this type of investigation, to develop methods of handling the great diversity of physical conditions encountered by cosmic rays propagating through high-energy emitting environments. Previously (Fatuzzo et al. 2010; Fatuzzo & Melia 2012b), we reported some results of this study pertaining to the spatial diffusion coefficients. We found that the spatial diffusion of particles through turbulent fields is not sensitive to the minimum wavelength of the fluctuations, so long as the particle's radius of gyration exceeds λ_min. For a given environment, as defined by B₀ and v_A, the diffusion process is thus dependent upon the maximum turbulence wavelength λ_max, the turbulent field strength, as characterized by η, and the turbulence spectrum, as characterized by the spectral index Γ.

We also noted that quasi-linear theory does not appear to be valid in the strong turbulence limit (see also O'Sullivan et al. 2009). We therefore investigated how the energy diffusion coefficient depends upon λ_max and η for Kolmogorov (Γ = 5/3) turbulence, and found that the energy diffusion coefficients could be characterized as $D_\gamma \propto \lambda _{{\rm max}}^{-0.47}$. This behavior is not consistent with quasi-linear theory, which instead predicts that $D_\gamma \propto \lambda _{{\rm max}}^{-0.67}$ for Kolmogorov turbulence in the strong turbulence (η ≳ 1) limit. However, we also found that D_γ∝η^1.2 in both the weak and strong turbulence limits.

Clearly, this initial sampling of the complex behavior of D_γ under a variety of physical conditions is far from satisfactory. The purpose of the present paper has been to complete this work, finding scaling relations that one may use to calculate D_γ under most conditions of interest, for all practical ranges of magnetic intensity, turbulence strength, ambient particle density, fluctuation size, and turbulence spectrum.

The empirical relations useful for this purpose have all been presented in Section 5. Broadly speaking, we have found that insofar as the energy diffusion coefficient D_γ is concerned, quasi-linear theory predicts its correct energy dependence only for very weak, isotropic turbulence (i.e., η ≲ 0.01). These predictions deviate substantially for η ∼ 1, particularly for turbulence spectral indices Γ > 1, such as Kolmogorov turbulence, which seems to be prevalent in many diverse environments. For example, for the physical conditions one encounters at the galactic center (i.e., B ∼ 1–10³ μG and n ∼ 1–10³ cm⁻³), the actual index characterizing the dependence of D_γ on γ may be as small as ∼1.4 instead of the predicted value ∼1.7.

In addition, our results indicate that there is no difference in the energy diffusion of particles between isotropic turbulence and the anisotropic turbulence profiles predicted by Goldreich & Sridhar (1995) so long as the turbulence is strong (δB ∼ B). On the other hand, our results indicate that anisotropic weak turbulence is considerably less effective in energetically scattering particles, consistent with the results of Chandran (2000).

Needless to say, if one is attempting to interpret the GeV Fermi spectrum of, say, the Fermi bubbles, in terms of an underlying population of cosmic rays, deviations from quasi-linear theory are rather critical, since the inferred particle distribution differs considerably from its injection point to where it emits the radiation, and the difference will be interpreted incorrectly with the inaccurate energy dependence predicted by these other techniques.

Not being sure of the normalization of D_γ has its own challenges. For one thing, it is not possible to say anything definitive about the overall power being generated by the acceleration of these cosmic rays, which speaks directly to the mechanism associated with the relativistic particle injection, or even to the required density of dark matter particles, if these cosmic rays are produced via dark matter decays and collisions. With our approach, it is not necessary to estimate the normalization of D_γ, because its absolute value is determined self-consistently from the statistical aggregate of numerous individual particle trajectories.

We have found that quasi-linear theory provides an acceptable estimate of the normalization of D_γ at ∼1 TeV energies, but can deviate considerably at lower energies, especially in the GeV range, and at energies exceeding 10–100 TeV. A large factor responsible for these differences is the incorrect energy dependence predicted for D_γ. Obviously, if the normalization is adequate at ∼1 TeV, the incorrect energy index will cause deviations at lower and higher energies.

Most importantly, however, our analysis has provided a method of determining not only the dependence of D_γ on γ, but also its absolute value without the need to normalize it from the data. Having said this, one is not completely free of ambiguity, since one must still have an accurate estimate of the physical conditions, i.e., the magnetic field, the ambient density, and other characteristics that determine the state of the medium though which the cosmic rays are propagating. Fortunately, these conditions are easier to measure than the diffusion coefficients themselves, and in a more sophisticated use of our technique, in which MHD turbulence is simulated numerically rather than via the simple Kolmogorov or Bohm scaling relations, one can approach a level of realism not available to any of the other methods.

We look forward to the application of the scaling relations we have presented in this paper across a broad range of physical environments, allowing us to study high-energy sources with a level of accuracy commensurate with the detailed measurements now being made by the ever improving suite of space-based and ground-based observatories.

This work was supported by Xavier University through the Hauck Foundation, and by ONR grant N00014-09-C-0032 at the University of Arizona.