Solutions of the cosmic ray velocity diffusion equation

Accepted Manuscript Solutions of the Cosmic Ray Velocity Diffusion Equation J. Lasuik, A. Shalchi PII: DOI: Reference:

S0273-1177(17)30461-1 http://dx.doi.org/10.1016/j.asr.2017.06.035 JASR 13289

To appear in:

Advances in Space Research

Received Date: Accepted Date:

1 May 2017 19 June 2017

Please cite this article as: Lasuik, J., Shalchi, A., Solutions of the Cosmic Ray Velocity Diffusion Equation, Advances in Space Research (2017), doi: http://dx.doi.org/10.1016/j.asr.2017.06.035

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Solutions of the Cosmic Ray Velocity Diffusion Equation J. Lasuik and A. Shalchi Department of Physics and Astronomy, University of Manitoba, Winnipeg, Manitoba R3T 2N2, Canada

Abstract In order to describe the propagation and acceleration of cosmic rays, one usually uses a diffusive transport equation. The most fundamental equation is the pitch-angle dependent diffusion equation which is usually called the Fokker-Planck equation. In the current paper we solve the position integrated equation numerically and analytically. For a constant pitch-angle Fokker-Planck coefficient we derive an exact solution of the corresponding transport equation and compare it with numerical solutions. We show that even if the scattering coefficient is assumed to be constant, the solution behaves well. As a second example we consider the case of a linear pitch-angle Fokker-Planck coefficient. Again we solve the corresponding transport equation numerically and analytically. In all cases considered, we find similar distribution functions. We also compute the corresponding velocity correlation functions and parallel diffusion coefficients. Our results are relevant for improving analytical theories of perpendicular diffusion, for code tests, and for different astrophysical applications where a pitch-angle dependent description of the particle motion is required. Key words: magnetic fields, turbulence, energetic particles

1

Introduction

A fundamental problem in astrophysics is to understand the origin and dynamics of cosmic rays. Due to complicated interactions between such energetic particles and astrophysical plasmas, their motion is usually of diffusive nature. Therefore, the propagation and acceleration of cosmic rays is described by a diffusive transport equation (see, e.g,. Schlickeiser 2002 and Zank 2014). It has Email address: [email protected] (J. Lasuik and A. Shalchi).

Preprint submitted to Elsevier

26 June 2017

to be emphasized, however, that in more recent years non-diffusive transport was discussed extensively in the literature (see, e.g., Zimbardo et al. 2006, Pommois et al. 2007, and Zimbardo et al. 2012). Due to the complexity of cosmic ray transport equations, full analytical solutions are usually not available and one has to perform numerical calculations. Usual transport equations are derived from the pitch-angle dependent cosmic ray Fokker-Planck equation (see, e.g., Schlickeiser 2002). Therefore, the Fokker-Planck equation is more fundamental compared to the usual diffusive transport equation. It is the purpose of the current paper to derive solutions for different forms of the pitch-angle scattering coefficient Dμμ and to compare those with solutions derived previously. We also solve the corresponding transport equations numerically to test the validity of our solutions. Furthermore, we compute particle velocity correlation functions and parallel diffusion coefficients. The remainder of this article is organized as follows. In Sect. 2 we briefly discuss the two-dimensional cosmic ray Fokker-Planck equation. In Sect. 3 we summarize solutions derived before in the literature. In Sects. 4 and 5 we solve the position integrated transport equation for a constant and a linear scattering coefficient, respectively. In Sect. 6 we present numerical solutions to the Fokker-Planck equations and compare them with our analytical findings. In Sect. 7 we summarize and conclude.

2

The Fokker-Planck Equation

The two-dimensional Fokker-Planck equation of cosmic ray transport is given by (see, e.g., Schlickeiser 2002 and Shalchi 2009) 

∂f ∂f ∂ ∂f + vμ = Dμμ (μ) ∂t ∂z ∂μ ∂μ



(1)

where we have used the time t, the position of the particle along the mean magnetic field z, the constant particle speed v, as well as the pitch-angle cosine μ = vz /v. The term on the right hand side of Eq. (1) describes pitchangle scattering and the scattering parameter Dμμ is usually called the pitchangle Fokker-Planck coefficient. In general, cosmic ray transport equations contain additional terms such as perpendicular diffusion, adiabatic focusing, or stochastic acceleration (see, e.g., Schlickeiser 2002). For some applications 2

one could just be interested in the z-integrated distribution function +∞ 

dz f (z, μ, t)

F (μ, t) =

(2)

−∞

which indicates the probability that the particle has the pitch-angle cosine μ at time t. By integrating Eq. (1) over z, we can easily derive 



∂ ∂F ∂F = Dμμ (μ) . ∂t ∂μ ∂μ

(3)

The form of F (μ, t) is entirely controlled by the pitch-angle Fokker-Planck coefficient Dμμ (μ). Any physical solution of Eq. (3) must satisfy the normalization condition 1 2

+1

dμ F (μ, t) = 1.

(4)

−1

For the initial distribution we assume F (μ, t = 0) = 2δ (μ − μ0 )

(5)

simply meaning that the particle has the initial pitch-angle cosine μ0 . It was shown in several papers and books, that if the Fokker-Planck equation (1) is averaged over the pitch-angle and a late time limit is considered, one obtains a diffusion equation (see, e.g., Schlickeiser 2002). The parallel diffusion coefficient κ therein is then related to the scattering coefficient Dμμ (μ) via the famous relation (see Earl 1974) v2 κ = 8

+1

dμ −1

(1 − μ2 )2 . Dμμ (μ)

(6)

In order to obtain a finite diffusion coefficient the integral in (6) has to be convergent. Furthermore, in order to derive the diffusion equation from (1), the following condition is usually assumed to be valid Dμμ (μ = ±1) = 0.

(7)

For most analytical forms proposed for the scattering parameter Dμμ , this condition is satisfied (see, e.g., Schlickeiser 2002 and Shalchi et al. 2009). 3

However, relation (6) can also be derived by using the less strict condition (see Shalchi 2011) 

∂F Dμμ (μ) ∂μ



= 0.

(8)

μ=±1

This condition is used in the current paper for investigations based on a constant scattering coefficient.

3

Previously Discussed Solutions

In the following we present different forms of Dμμ discussed previously in the literature. In Sect. 6 the corresponding solutions of Eq. (3) will be compared with the new solutions derived in the current paper. In Fig. 1 all models used in the current paper for Dμμ are visualized.

3.1 The Scatter-Free Solution

As an illustrative example we consider the simplest scenario, namely the scatter-free case corresponding to Dμμ (μ) = 0. For this case, Eq. (3) becomes ∂F =0 ∂t

(9)

and the solution is any time-independent function of μ. Since we need to satisfy the initial condition (5), the solution for the scatter-free case is, of course, F (μ, t) = 2δ (μ − μ0 )

(10)

meaning that the particle has always the same pitch-angle.

3.2 An Isotropic Scattering Coefficient

A simple model for the scattering coefficient Dμμ is the isotropic model in which we have by definition 



Dμμ (μ) = D 1 − μ2 .

(11)

4

It was shown that the isotropic form should be valid in cases where the turbulent magnetic field is stronger than the mean field (see Shalchi et al. 2009). Previously the averaged Fokker-Planck equation (3) was solved for this specific form of the scattering coefficient (see Shalchi 2006). The following solution was found F (μ, t) = 1 +

∞ 

(2l + 1)Pl (μ0 )Pl (μ)e−l(l+1)Dt

(12)

l=1

where we have used the Legendre polynomials Pl (μ). For t → ∞, we find F (μ, t → ∞) → 1 corresponding to a pitch-angle isotropization process. For t → 0, on the other hand, Eq. (12) becomes equivalent to Eq. (10) as required. The parallel velocity correlation function of the energetic particles is defined via v2 vz (t)vz (0) = 4

+1

+1

dμ −1

dμ0 μμ0 F (μ, t) .

(13)

−1

By using Eq. (12) for the function F (μ, t) and the orthogonality relation of the Legendre polynomials +1

dx Pn (x)Pm (x) = −1

2 δmn 2m + 1

(14)

we can easily derive vz (t)vz (0) =

v 2 −2Dt e . 3

(15)

A diffusion coefficient can be calculated by employing the Taylor-Green-Kubo formula (see Taylor 1922, Green 1951, and Kubo 1957) ∞

κ =

dt vz (t)vz (0) .

(16)

0

It was shown in Shalchi (2011) that Eq. (6) can be derived by combining the Fokker-Planck equation (1) with the latter formula. This derivation can be performed for an arbitrary pitch-angle FokkerPlanck coefficient Dμμ . 5

For the isotropic form (15), the Taylor-Green-Kubo formula yields κ =

v2 . 6D

(17)

We can derive the same formula by combining Eq. (11) with (6). For the case considered here, the parallel mean free path is given by 3 v λ = κ = . v 2D

(18)

Therefore, the velocity correlation function (15) can be written as vz (t)vz (0) =

v 2 −vt/λ e . 3

(19)

The latter form is often used in non-linear theories for perpendicular diffusion in order to model the parallel motion of the particles (see, e.g., Owens 1974 and Matthaeus et al. 2003). However, this simple exponential form is only exact as long as pitch-angle scattering is isotropic (see Shalchi 2011).

3.3 The approximated quasi-linear case

Previously, the quasi-linear Fokker-Planck coefficient was discussed for the case that the inertial range spectral index is s = 2 (see Shalchi 2011). In this case we have Dμμ = D |μ| (1 − μ2 ) ≡ D |μ| (1 − |μ|)(1 + |μ|).

(20)

Since 1 + |μ| is just a number between 1 and 2, this form can be approximated by Dμμ ≈ D |μ| (1 − |μ|) .

(21)

This form has the same important properties as the quasi-linear form (20), e.g., Dμμ (μ = ±1) = 0 and Dμμ (μ = 0) = 0. The latter condition can cause problems if the quasi-linear Fokker-Planck coefficient is used in Eq. (6). This problem, and the fact that quasi-linear theory does in general not work for pitch-angles at and close to 90o , is well-known in diffusion theory as the 90o problem. A solution to this problem was proposed in Shalchi (2005) and Tautz et al. (2008). 6

1.2

1

Dμμ / D

0.8

0.6

0.4

0.2

0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

μ

Fig. 1. The different models considered in the current paper for the pitch-angle scattering coefficient Dμμ . Shown are the constant coefficient (dotted line), the linear model (dashed line), the approximated quasi-linear model (dash-dotted line), and the isotropic model (solid line).

The averaged equation (3) has been solved for the approximated quasi-linear case in Shalchi (2011) and the following solution was derived

F ± (μ, t) = 2Θ (±μ0 )

∞ 

(2n + 1)Pn (x0 = 1 ∓ 2μ0 )

n=0

× Pn (x = 1 ∓ 2μ)e−n(n+1)Dt .

(22)

The solution F + (μ, t) has to be used for positive pitch-angle cosine values and F − (μ, t) for negative values. Here we have used the Heaviside step function Θ(x). The meaning of this function in Eq. (22) is the following. If μ0 < 0, the function F + is zero due to the step function and only F − is not zero. For μ0 > 0, on the other hand, we have F + = 0 but F − = 0. This means that the particle cannot be scattered through μ = 0 because of Dμμ (μ = 0) = 0. As shown before the velocity correlation function for the scattering coefficient used here is (see Shalchi 2011) 



v2 1 1 + e−2Dt . vz (t)vz (0) = 4 3

(23)

After employing the Taylor-Green-Kubo formula (16), we can clearly see that we obtain super-diffusive transport. Furthermore, Eq. (6) does not convergence in this case. 7

4

A Constant Scattering Coefficient

The simplest non-trivial case corresponds to a constant pitch-angle FokkerPlanck coefficient Dμμ (μ) = D = const.

(24)

In this case the integral (6) is convergent but this form of Dμμ does not satisfy condition (7). However, we can chose the solution so that condition (8) is satisfied. In the following we solve Eq. (3) for a constant Dμμ and thereafter we compute velocity correlation functions and the parallel diffusion coefficient.

4.1 The solution F (μ, t) For a constant scattering coefficient, Eq. (3) becomes ∂2F ∂F =D 2. ∂t ∂μ

(25)

The latter equation is a typical diffusion equation which can be solved by using standard methods. However, one has to be careful because the variable μ is limited due to −1 ≤ μ ≤ +1. As a first step we employ the Ansatz F (μ, t) = P (μ) ρ (t)

(26)

to find 1 ∂2P 1 ∂ρ =D . ρ ∂t P ∂μ2

(27)

To proceed we introduce the constant ω so that we derive the two ordinary differential equations ∂ρ = −ω 2 Dρ ∂t

(28)

∂2P = −ω 2 P. ∂μ2

(29)

and

8

The two differential equations can easily be solved by ρ (ω, t) ∝ e−ω

2 Dt

(30)

and P (ω, μ) = a (ω) cos (ωμ) + b (ω) sin (ωμ) .

(31)

Therefore, the general solution of Eq. (25) is given by F (μ, t) =



2

Pn (μ) e−ωn Dt .

(32)

n

The latter solution must satisfy the normalization condition (4). In the case considered here, Eq. (8) becomes 

∂ F (μ, t) ∂μ



= 0.

(33)

μ=±1

One can easily show that if the latter conditions are combined with Eq. (25), the normalization condition (4) is automatically satisfied. Because of Eq. (33), we find 



∂ Pn (μ) ∂μ

μ=±1

= [−a (ω) ω sin (ωμ) + b (ω) ω cos (ωμ)]μ=±1 = 0.

(34)

Therefore, we can find the two conditions b (ω) cos (ω) + a (ω) sin (ω) = 0

(35)

b (ω) cos (ω) − a (ω) sin (ω) = 0.

(36)

and

By adding and subtracting the latter two equations, this becomes b (ω) cos (ω) = 0 a (ω) sin (ω) = 0.

(37) 9

For a = 0 and b = 0 this means that ω = π(n + 1/2) and for a = 0 and b = 0 we need to have ω = πn where n = 0, 1, 2, . . .. Therefore, the solution of our problem can be written as

F (μ, t) = +

∞  n=0 ∞ 

2 Dt

an cos (πnμ) e−(πn)

2 Dt

bn sin (π(n + 1/2)μ) e−(π(n+1/2))

.

n=0

(38) To determine the constants an and bn , the following relations are useful (see Jeffrey & Dai 2008) +1

dμ sin (nπμ) sin (mπμ) = δmn −1 +1

dμ cos (nπμ) cos (mπμ) = δmn

(39)

−1

and +1

dμ sin (aμ) cos (bμ) = 0.

(40)

−1

For t = 0, Eq. (38) becomes

F (μ, t = 0) = +

∞  n=0 ∞ 

an cos (πnμ) bn sin (π(n + 1/2)μ)

n=0

= 2δ (μ − μ0 )

(41)

where we have used the initial condition (5). Now we multiply the latter equation by cos(πmμ) and integrate over all μ to find +1

2

dμ cos (πmμ) δ (μ − μ0 )

−1

=

∞  n=0

+1

an

dμ cos (πnμ) cos (πmμ) −1

10

+

∞ 

+1

bn

n=0

dμ sin (π(n + 1/2)μ) cos (πmμ) .

(42)

−1

By employing Eqs. (39) and (40), one can easily show that the second integral is zero and, thus, we derive am = 2 cos (πmμ0 ) .

(43)

Now we multiply Eq. (41) by sin(π(m + 1/2)μ) and integrate over all μ +1

2

dμ sin (π(m + 1/2)μ) δ (μ − μ0 )

−1

=

+

∞ 

+1

an

dμ cos (πnμ) sin (π(m + 1/2)μ)

n=0

−1

∞ 

+1

n=0

bn

dμ sin (π(n + 1/2)μ) sin (π(m + 1/2)μ) . −1

(44) By employing again Eqs. (39) and (40), one can easily show that the first integral is zero and, thus, we derive bm = 2 sin (π(m + 1/2)μ0 ) .

(45)

Finally the solution of Eq. (25) can be written as

F (μ, t) = 1 + 2 +2

∞ 

∞ 

2 Dt

cos (πnμ0 ) cos (πnμ) e−(πn)

n=1

sin (π(n + 1/2)μ0 )

n=0 2 Dt

× sin (π(n + 1/2)μ) e−(π(n+1/2))

.

(46)

We like to point out that our solution can be expressed by Jacobi theta functions but this is not useful for the investigations present in the current paper. We can recover the scatter-free solution by setting D = 0 in Eq. (46). This gives us F (μ, t = 0) 11

=1 + 2 +2

∞ 

∞ 

cos (πnμ0 ) cos (πnμ)

n=1

sin (π(n + 1/2)μ0) sin (π(n + 1/2)μ)

n=0

= 2δ (μ − μ0 )

(47)

as expected. From Eq. (46) we can easily derive F (μ, t → ∞) → 1 corresponds to a pitch-angle isotropization process.

4.2 Velocity correlation function

The velocity correlation function can be calculated via Eq. (13). Together with solution (46) we find vz (t)vz (0) ∞ +1 v2  = dμ0 μ0 cos (πnμ0 ) 2 n=1 −1

+1

2 Dt

dμ μ cos (πnμ) e−(πn)

×

−1 ∞ +1 v2  + dμ0 μ0 sin (π(n + 1/2)μ0 ) 2 n=0 −1

+1

2 Dt

dμ μ sin (π(n + 1/2)μ) e−(π(n+1/2))

×

.

(48)

−1

The occurring integrals are solved by +1

dμ μ cos (πnμ) = 0

(49)

−1

and

1

+1

dμ μ sin (π(n + 1/2)μ) = −1

8 (−1)n . π 2 (1 + 2n)2

(50)

+1 d The integral can be written as − dα −1 dμ cos(αμ). The integral over the cosine function is straightforward to evaluate. After calculating the derivative with respect to α and by setting α = π(n + 1/2) we find the solution of the integral. 1

12

Therewith, the velocity correlation function can be written as vz (t)vz (0) =

∞ 1 32v 2  −(π(n+1/2))2 Dt . 4e 4 π n=0 (1 + 2n)

(51)

In the late time limit the main contribution comes from n = 0 and, therefore, vz (t → ∞)vz (0) →

32v 2 −(π/2)2 Dt e . π4

(52)

For t = 0, on the other hand, we have 



vz2 (0) =

∞ 1 32v 2  . 4 π n=0 (1 + 2n)4

(53)

The sum yields π 4 /96 (see Jeffrey & Dai 2008) and, therefore, we find 



vz2 (0) =

v2 . 3

(54)

We can clearly see that the velocity correlation function obtained for a constant scattering coefficient (51) can be written as a sum over exponential functions. The simple form (15) obtained for an isotropic pitch-angle FokkerPlanck coefficient can no longer be found. However, we like to point out that Eq. (15) would still be a good approximation for Eq. (51). In Fig. 10 we have shown the velocity correlation function obtained here for a constant Dμμ in comparison with the results obtained above for an isotropic scattering coefficient (15).

4.3 Parallel diffusion coefficient

A diffusion coefficient can be calculated by employing the Taylor-Green-Kubo formula (see Eq. (16) of the current paper). For the velocity correlation function we employ Eq. (51) to find ∞ ∞ 32v 2  1 2 κ = 4 dt e−((n+1/2)π) Dt 4 π n=0 (1 + 2n) 0

=

∞ 128v 2  1 6 π D n=0 (1 + 2n)6

(55)

13

The sum yields π 6 /960 (see Jeffrey & Dai 2008) and, therefore, we find for the parallel diffusion coefficient κ =

2v 2 15D

(56)

and for the parallel mean free path λ =

2v . 5D

(57)

We like to point out, that these formulas agree perfectly with Eq. (6). This is an important result because when Eq. (6) is derived, it is assumed that Dμμ (μ = ±1) = 0. This condition is not satisfied for a constant Dμμ = D. We have shown, however, that even in this case, Eq. (6) provides the correct result. Most probably this is the case because we still satisfy condition (8).

5

A Linear Scattering Coefficient

The true form of the scattering parameter Dμμ is not known and can be very complicated (see, e.g., Qin & Shalchi 2009 for an investigation of this matter based on test-particle simulations). For applications such as solving the Fokker-Planck equation one desires simple analytical forms. The constant form used above, however, could be an oversimplification because it is known that the real scattering coefficient does depend on the pitch-angle cosine μ. In the current section we employ an alternative form to explore the importance of the pitch-angle dependence and its influence on the solution of the corresponding Fokker-Planck equation. This alternative form is linear meaning that we assume Dμμ (μ) = D (1 − |μ|) .

(58)

Compared to the constant coefficient used above, this linear form has the advantage, that condition (7) is now satisfied. In the following we solve Eq. (3) for this case and compute velocity correlation functions and the parallel diffusion coefficient. 14

5.1 The solution F (μ, t) Here we have to distinguish between positive and negative μ. For the former case, the position integrated Fokker-Planck equation (3) becomes 



∂ ∂F + ∂F + =D (1 − μ) . ∂t ∂μ ∂μ

(59)

Here F + denotes that the solution can only be used for positive values of μ. As a first step we employ the Ansatz F + (μ, t) = P + (μ) ρ (t)

(60)

to find 



∂P + 1 ∂ρ 1 ∂ (1 − μ) . =D + ρ ∂t P ∂μ ∂μ

(61)

To proceed we introduce the constant ω so that we derive the two ordinary differential equations ∂ρ = −ω 2 Dρ ∂t

(62)

and 



∂P + ∂ (1 − μ) = −ω 2 P + . ∂μ ∂μ

(63)

The first differential equation can easily be solved by ρ (ω, t) ∝ e−ω

2 Dt

.

(64)

The second equation, however, is more difficult to solve. First we employ the transformation

x = 2ω 1 − μ.

(65)

After lengthy straightforward algebra, we obtain x2

∂P + ∂2P + + x2 P + = 0 + x ∂x2 ∂x

(66)

15

corresponding to Bessel’s differential equation. Solutions are, of course, first and second kind Bessel functions. The second kind or Weber functions Yν are not part of our solution because they have a singularity at the origin x = 0. Therefore, our solution is given by the function J0 (x). Therefore, the general solution of Eq. (59) has the form +

F (μ, t) =

 n

αn+ J0





2

2ωn 1 − μ e−ωn Dt .

(67)

Next we consider the case of negative values of μ. We can repeat all the steps performed above in order to solve 



∂F − ∂F − ∂ (1 + μ) . =D ∂t ∂μ ∂μ

(68)

Steps (60)-(62) are as above but we find the following equation instead of (63) 



∂P − ∂ (1 + μ) = −ω 2 P − . ∂μ ∂μ

(69)

Now we employ the transformation

x = 2ω 1 + μ

(70)

leading again to the Bessel differential equation (66). Therefore, the solution for negative μ has the form F − (μ, t) =

 n





2

αn− J0 2ωn 1 + μ e−ωn Dt .

(71)

Our solution must satisfy the normalization condition (4). By employing Eqs. (67) and (71) we can write

1= ×

 1  + αn + αn− 2 n 1





2

dμ J0 2ωn 1 − μ e−ωn Dt .

(72)

0

In order to rewrite the μ-integral, we use transformation

x=

1−μ

(73)

16

to find 1





dμ J0 2ω 1 − μ = 2

0

1

dx xJ0 (2ωx) 0

1 = J1 (2ω) . ω

(74)

Therewith, condition (72) becomes

1=

 1 1  + 2 αn + αn− J1 (2ωn ) e−ωn Dt . 2 n ωn

(75)

Since the left hand side of this equation is constant, the right hand side must be time-independent as well. This is only possible if ω = 0 and if the right hand side is zero for ω = 0. For the former case the corresponding coefficients are α0+ = α0− = 1.

(76)

One possibility to have the right hand side equal to zero is that the Bessel function itself is zero. We define the zeros zn of the corresponding Bessel function via J1 (zn ) = 0.

(77)

Therefore, the Eigenvalues for this case are 1 ωn = zn . 2

(78)

In Table 1 we list the first few values of the zeros zn . We also need to satisfy the connection condition F + (μ = 0, t) = F − (μ = 0, t)

(79)

and, therefore, we need for this case αn+ = αn− =: αnS

(80)

where the index S stands for symmetric because this part of the solution is symmetric in μ. Another way of satisfying condition (75) is αn+ = −αn− =: αnA

(81) 17

where the index A stands for anti-symmetric because this part of the solution is anti-symmetric in μ. We also have to satisfy condition (79). Therefore, the corresponding Eigenvalues need to satisfy 1 ωn = y n 2

(82)

where we have used the zeros of the Bessel function defined via J0 (yn ) = 0.

(83)

In Table 1 we also list the first few values of the zeros yn . Therefore, our solutions have the form

+

F (μ, t) = 1 + +

∞ 

∞  n=1



αnS J0 



1 2

zn 1 − |μ| e− 4 zn Dt

1 2

αnA J0 yn 1 − |μ| e− 4 yn Dt

(84)

n=1

for positive μ and

−

F (μ, t) = 1 + −

∞ 

∞  n=1



αnS J0 



1 2

zn 1 − |μ| e− 4 zn Dt



1 2

αnA J0 yn 1 − |μ| e− 4 yn Dt

(85)

n=1

for negative μ. In order to determine the remaining coefficients αnS and αnA , we use the initial distribution (5).

5.1.1 The coefficients αnS In this case we have to use solution (84) for the initial time t = 0. Then we multiply the latter condition with a Bessel function with a different n and integrate to find +1





dμ δ (μ − μ0 ) J0 zm 1 − |μ|

−1

18

Table 1 The first 20 zeros of the Bessel functions used in the current paper. We define the zeros via J0 (yn ) = 0 and J1 (zn ) = 0. The values were taken from the literature (see Abramowitz & Stegun 1974).

=

∞ 

1

αnS

n=1

=2

∞  n=1



n

yn

zn

1

2.4048255577

3.83171

2

5.5200781103

7.01559

3

8.6537279129

10.17347

4

11.7915344391

13.32369

5

14.9309177086

16.47063

6

18.0710639679

19.61586

7

21.2116366299

22.76008

8

24.3524715308

25.90367

9

27.4934791320

29.04683

10

30.6346064684

32.18968

11

33.7758202136

35.33231

12

36.9170983536

38.47477

13

40.0584257646

41.61709

14

43.1997917132

44.75932

15

46.3411883717

47.90146

16

49.4826098974

51.04354

17

52.6240518411

54.18555

18

55.7655107550

57.32753

19

58.9069839261

60.46946

20

62.0484691902

63.61136





dμ J0 zn 1 − μ J0 zm 1 − μ



0

1

αnS

dx xJ0 (zn x) J0 (zm x) 0

(86) where the anti-symmetric part of the solution automatically disappeared. Furthermore we have used again transformation (73). The first integral is trivial due to the Dirac delta. In order to solve the x-integral, we employ (see 19

Abramowitz & Stegun 1974) 1 0

1 dx xJ0 (zn x) J0 (zm x) = δmn J02 (zm ) . 2

(87)

Therewith Eq. (86) becomes 



S 2 J0 (zm ) J0 zm 1 − |μ0 | = αm

(88)

and we find for the coefficients



αnS =

J0 zn 1 − |μ0 |



J02 (zn )

.

(89)

Therefore, the coefficients in the symmetric part of the solution are finally known.

5.1.2 The coefficients αnA To determine the remaining coefficients, we consider the following quantity

1 2 −

=

1 2

1

+







dμ F (μ, t = 0) J0 ym 1 − |μ|

0

0

dμ F − (μ, t = 0) J0 ym 1 − |μ|

−1

∞ 

1

αnA

n=1

=2



∞  n=1







dμ J0 yn 1 − |μ| J0 ym 1 − μ



0

1

αnA

dx xJ0 (yn x) J0 (ym x)

(90)

0

where we have used Eqs. (84) and (85) and we have employed again transformation (73). To solve the x-integral, we use 1 0

1 dx xJ0 (yn x) J0 (ym x) = δmn J12 (ym ) . 2

20

(91)

For the first integrals in Eq. (90) we employ again initial condition (5). The first integral gives only a contribution if μ0 ≥ 0 and the second one only if μ0 ≤ 0 due to the Dirac delta. Therefore, we find



αnA (μ0 ≥ 0) =



J0 yn 1 − |μ0 |

(92)

J12 (yn )

and 

αnA (μ0 ≤ 0) = −

J0 yn 1 − |μ0 |



(93)

J12 (yn )

and the coefficients αnA are known. 5.1.3 Final solution of the problem If we combine the coefficients (89), (92), and (93) with solutions (84) and (85) we finally obtain F (μ, t) =1 + +



∞ 

J0 (zn 1 − |μ0 |)J0 (zn 1 − |μ|)

n=1

J02 (zn )

∞ 

μ0 μ |μ0 | |μ| n=1

1 2

e− 4 zn Dt

J0 (yn 1 − |μ0 |)J0 (yn 1 − |μ|) J12 (yn )

1 2

e− 4 yn Dt (94)

for arbitrary −1 ≤ μ ≤ 1. The fraction in front of the second sum is just there to ensure that we have the correct sign. In the following we use this solution to determine velocity correlation functions and the parallel spatial diffusion coefficient as it was done above for the constant scattering parameter.

5.2 Velocity correlation function

By combining Eq. (13) with solution (94), we derive for the velocity correlation function vz (t)vz (0) = v 2

∞ 

1

2 n=1 J1 (yn )

1 2

e− 4 yn Dt

21

⎡ 1  × ⎣ dμ



⎤ 2 1−μ ⎦ .

μJ0 yn

(95)

0

The integral therein can be solved by 1





dμ μJ0 yn 1 − μ =

0

4 J2 (yn ) yn2

(96)

and the velocity correlation function can be written as vz (t)vz (0) = 16v 2

∞ 

J22 (yn ) − 1 yn2 Dt e 4 . 4 J 2 (y ) y n n=1 n 1

(97)

Now we employ the relation (see, e.g., Abramowitz & Stegun 1974) Jν (z) =

z [Jν−1 (z) + Jν+1 (z)] 2ν

(98)

which becomes in our case J1 (yn ) =

yn J2 (yn ). 2

(99)

Therefore, we finally obtain vz (t)vz (0) = 64v

2

∞ 

1 − 1 yn2 Dt e 4 . 6 y n n=1

(100)

Again the velocity correlation function is a sum over exponential functions and the simple form (19) often used in diffusion theory is no longer valid. For t = 0 this becomes 

2

vz (0)



= 64v

2

∞ 

1 . 6 n=1 yn

(101)

According to Sneddon (1960) we have ∞ 

1 1 = 6 192 n=1 yn

(102)

22

and, thus, 



vz (0)2 =

v2 3

(103)

which is the result we have expected. For late times, only the first term in Eq. (100) contributes to the sum and we can approximate vz (t → ∞)vz (0) →

64 2 − 1 y12 Dt v e 4 . y16

(104)

Again we find an exponential function for the velocity correlation function in the late time limit. 5.3 Parallel diffusion coefficient

In order to compute the parallel diffusion coefficient we combine the TaylorGreen-Kubo formula (16) with Eq. (100) to derive ∞ 1  2 Dt − 14 yn κ = 64v dt e 6 n=1 yn 2

=

256v D

∞ 

∞ 2 

0

1 . 8 n=1 yn

(105)

To continue we use (see Sneddon 1960) ∞ 

1 11 . = 8 12288 n=1 yn

(106)

Therefore, we finally find for the parallel spatial diffusion coefficient 11v 2 κ = 48D

(107)

and for the parallel mean free path λ =

11v . 16D

(108)

Based on Eq. (6), we find exactly the same results. 23

6

Numerical results and visualization

In Sects. 3-5 of the current paper we presented analytical solutions to the Fokker-Planck equation (3). Alternatively, we can solve this type of differential equation numerically. In the following we solve the more general transport equation (1) numerically and then we use the obtained function f (z, μ, t) to compute the z-integrated function F (μ, t) as well as the running diffusion coefficient d(t) defined via d(t) =

1  (Δz)2 . 2t

(109)

All our numerical findings are compared with the analytical results derive above. In Figs. 2-5 we show the numerically obtained function F (μ, t) for different times. For the initial pitch-angle cosine we have set μ0 = 0.5 as an example. For the constant, the linear, and the isotropic Dμμ we find very similar solutions. Very clearly we can observe a pitch-angle isotropization process. The approximated quasi-linear case is an exception. The solution behaves very differently because in this case particles cannot be scattered through μ = 0 due to Dμμ (μ = 0) = 0. In Figs. 6-9 we compare the numerical solutions with our analytical solutions for the constant and the linear Dμμ , respectively. Very clearly we can see that the agreement between analytical and numerical solutions is almost perfect. A small discrepancy can be observed for early times (see Fig. 9). This is because we used a (very narrow) Gaussian distribution for the initial function in our numerical work. In the analytical treatment an (exact) Dirac delta has been used (see Eq. (5) of the current paper). In Fig. 10 we have visualized the velocity correlation functions as given by Eqs. (19), (51), and (100). We plotted the corresponding functions versus the dimensionless time vt/λ . According to Fig. 10, all velocity correlation functions are very similar confirming the validity of model (19) often used in non-linear theories of perpendicular diffusion. In Fig. 11 we have shown the running diffusion coefficient as given by Eq. (109) versus time. Shown are the numerical findings for the constant and linear Dμμ , respectively. For the sake of comparison we have also shown the analytical results as given by Eqs. (56) and (107). We can clearly see that the numerical solution initially obeys a complicated time-dependence and then turns into a diffusive regime in perfect agreement with the analytical result. 24

8 7 6

F (μ)

5 4 3 2 1 0 -1

-0.5

0

μ

0.5

1

Fig. 2. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.01. Shown are the results for the isotropic Dμμ (dotted line), the constant Dμμ (dashed line), the approximated quasi-linear Dμμ (dash-dotted line), and the linear Dμμ (solid line).

7

Summary and Conclusion

In the current paper we have explored Eq. (3) which is the position integrated cosmic ray Fokker-Planck equation. The solution of this equation F (μ, t) describes the probability that the particle has the pitch-angle cosine μ at time t. In order to solve partial differential equation (3) one has to specify the scattering coefficient Dμμ also known as pitch-angle Fokker-Planck coefficient. In previous work this differential equation has been solved for an isotropic scattering coefficient as well as for an approximated quasi-linear Fokker-Planck coefficient. In the current paper we have solved this equation for two more cases, namely for a constant scattering parameter Dμμ = D and for the linear form Dμμ = D(1 − |μ|). These forms are visualized together with previously used models in Fig. 1. There are important applications for the work presented in the current paper. Some examples are: (1) To derive solutions of the Fokker-Planck equation as done in the current paper helps us to understand the transport of energetic particles. For 25

3.5 3

F (μ)

2.5 2 1.5 1 0.5 0 -1

-0.5

0

μ

0.5

1

Fig. 3. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.1. Shown are the results for the isotropic Dμμ (dotted line), the constant Dμμ (dashed line), the approximated quasi-linear Dμμ (dash-dotted line), and the linear Dμμ (solid line).

instance, we have computed velocity correlation functions along the mean magnetic field for different models of the scattering parameter Dμμ . Such velocity correlation functions are required in some non-linear theories for the perpendicular diffusion coefficient (see Owens 1974 and Matthaeus et al. 2003). Furthermore, we can directly test Eq. (6) for the parallel diffusion coefficient. (2) Our work is useful for code test. One can be interested in the numerically solution of more general and realistic versions of the cosmic ray FokkerPlanck equation (see Schlickeiser 2002). Such general codes can easily be tested by comparing numerical solutions for special forms of the scattering coefficient Dμμ with the exact solutions derived in the current paper. (3) The analytical solutions presented and discussed in the current paper can be used as starting point for more general investigations. For instance, one could try to solve the position dependent equation (1) or an equation containing additional terms such as adiabatic focusing. The methods and solutions derived in the present paper will be useful in order to perform such more complicated tasks. (4) Analytical solution to the Fokker-Planck equation as derived in the current paper can be used to estimate mean free paths of impulsive solar 26

2.5

2

F (μ)

1.5

1

0.5

0 -1

-0.5

0

μ

0.5

1

Fig. 4. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.5. Shown are the results for the isotropic Dμμ (dotted line), the constant Dμμ (dashed line), the approximated quasi-linear Dμμ (dash-dotted line), and the linear Dμμ (solid line).

energetic particle events by fitting the anisotropy profiles obtained from spacecraft observations (see, e.g., He & Qin 2011 and He & Wan 2012). Due to the work presented in the current paper, there are now four cases for which an exact analytical solution can be derived. In Table 2 we summarize those cases and refer to the corresponding solution. It will be subject of future work to use those solutions as starting point in order to explore the full position dependent Fokker-Planck equation analytically.

Acknowledgements

A. Shalchi acknowledges support by the Natural Sciences and Engineering Research Council (NSERC) of Canada. 27

3

2.5

F (μ)

2

1.5

1

0.5

0 -1

-0.5

0

μ

0.5

1

Fig. 5. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 5. Shown are the results for the isotropic Dμμ (dotted line), the constant Dμμ (dashed line), the approximated quasi-linear Dμμ (dash-dotted line), and the linear Dμμ (solid line). Table 2 Different solutions of the cosmic ray Fokker-Planck equation (3). The parameter Dμμ is the pitch-angle Fokker-Planck coefficient and D is a constant. Dμμ (μ)

Exact solution

Type of solution

Parallel diffusion coefficient

0

Eq. (10)

Dirac delta

Unperturbed motion

D

Eq. (46)

Trigonometric functions

κ = 2v 2 /(15D)

D(1 − |μ|)

Eq. (94)

Bessel functions

κ = 11v 2 /(48D)

D(1 − μ2 )

Eq. (12)

Legendre polynomials

κ = v 2 /(6D)

D |μ| (1 − |μ|)

Eq. (22)

Shifted Legendre polynomials

Superdiffusive motion

References

[1] Abramowitz M., Stegun I.A., 1974. Handbook of Mathematical Functions. Dover Publications, New York. [2] Earl J.A., 1974. The diffusive idealization of charged-particle transport in random magnetic fields. Astrophys. J. 193, 231-242.

28

8 7 6

F (μ)

5 4 3 2 1 0 -1 -1

-0.8

-0.6

-0.4

-0.2

0

μ

0.2

0.4

0.6

0.8

1

Fig. 6. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.01. Shown are the numerical solution for constant Dμμ (dotted line), the analytical solution for constant Dμμ (dashed line), the numerical solution for a linear Dμμ (dash-dotted line), and the analytical solution for a linear Dμμ (solid line). [3] Green, M. S., 1951. Brownian Motion in a Gas of Noninteracting Molecules. J. Chem. Phys. 19, 1036-1046. [4] He, H.-Q., Qin, G., 2011. A Simple Analytical Method to Determine Solar Energetic Particles’ Mean Free Path. Astrophys. J. 730, 46. [5] He, H.-Q., Wan, W., 2012. A Direct Method to Determine the Parallel Mean Free Path of Solar Energetic Particles with Adiabatic Focusing. Astrophys. J. 747, 38. [6] Jeffrey, A., Dai, H.-H., 2008. Handbook of Mathematical Formulas and Integrals, Elsevier Academic Press, London, New York, San Diego. [7] Kubo, R., 1957. Statistical-Mechanical Processes. J. Phys. Soc. Jpn. 12, 570-586.

Theory

of

Irreversible

[8] Matthaeus, W.H., Qin, G., Bieber J.W., Zank, G.P., 2003. Nonlinear Collisionless Perpendicular Diffusion of Charged Particles. Astrophys. J. 590, L53-L56. [9] Owens, A.J., 1974. The Effects of Nonlinear Terms in Cosmic-Ray Diffusion Theory. Astrophys. J. 191, 235-244. [10] Pommois, P., Zimbardo, G., Veltri, P., 2007. Anomalous, non-Gaussian transport of charged particles in anisotropic magnetic turbulence. Phys. Plasmas 14, 012311.

29

3

2.5

F (μ)

2

1.5

1

0.5

0 -1

-0.8

-0.6

-0.4

-0.2

0

μ

0.2

0.4

0.6

0.8

1

Fig. 7. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.1. Shown are the numerical solution for constant Dμμ (dotted line), the analytical solution for constant Dμμ (dashed line), the numerical solution for a linear Dμμ (dash-dotted line), and the analytical solution for a linear Dμμ (solid line). [11] Qin, G., Shalchi, A., 2009. Pitch-Angle Diffusion Coefficients of Charged Particles from Computer Simulations. Astrophys. J. 707, 6166. [12] Schlickeiser, R., 2002. Cosmic Ray Astrophysics, Springer, Berlin. [13] Shalchi, A., 2005. Second-order quasilinear theory of cosmic ray transport. Phys. Plasmas 12, 052905. [14] Shalchi, A., 2006. Analytical investigation of the two-dimensional cosmic ray Fokker-Planck equation. Astron. Astrophys. 448, 809-816. [15] Shalchi, A., 2009. Nonlinear Cosmic Ray Diffusion Theories. Astrophysics and Space Science Library, vol. 362. Springer, Berlin. ˇ [16] Shalchi, A., Skoda, T., Tautz, R. C., Schlickeiser, R., 2009. Analytical description of nonlinear cosmic ray scattering: isotropic and quasilinear regimes of pitch-angle diffusion. Astron. Astrophys. 507, 589-597. [17] Shalchi, A., 2011. Velocity correlation functions of charged particles derived from the Fokker-Planck equation. Adv. Space Res. 47, 1147-1164. [18] Sneddon, I.N., 1960. On some infinite series involving the zeros of Bessel functions of the first kind. Proc. Glasgow Math. Assoc. 4, 144156. [19] Tautz, R. C., Shalchi, A., Schlickeiser, R., 2008. Solving the 90 Scattering Problem in Isotropic Turbulence. Astrophys. J. 685, L165.

30

3.5 3

F (μ)

2.5 2 1.5 1 0.5 0 -1

-0.8

-0.6

-0.4

-0.2

0

μ

0.2

0.4

0.6

0.8

1

Fig. 8. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 0.5. Shown are the numerical solution for constant Dμμ (dotted line), the analytical solution for constant Dμμ (dashed line), the numerical solution for a linear Dμμ (dash-dotted line), and the analytical solution for a linear Dμμ (solid line). [20] Taylor, G., 1922. Diffusion by Continuous Movements. Proc. London Math. Soc. 20, 196-212. [21] Zank, G.P., 2014. Transport Processes in Space Physics and Astrophysics. Lecture Notes in Physics 877, Springer, New York. [22] Zimbardo, G., Pommois, P., Veltri, P., 2006. Superdiffusive and Subdiffusive Transport of Energetic Particles in Solar Wind Anisotropic Magnetic Turbulence. Astrophys. J. 639, L91-L94. [23] Zimbardo, G., Perri, S., Pommois, P., Veltri, P., 2012. Anomalous particle transport in the heliosphere. Adv. Space Res. 49, 1633-1642.

31

1.5

1

0.5 -1

-0.5

0

0.5

1

Fig. 9. Numerical solutions to F (μ, t) plotted versus μ at time τ = Dt = 5. Shown are the numerical solution for constant Dμμ (dotted line), the analytical solution for constant Dμμ (dashed line), the numerical solution for a linear Dμμ (dash-dotted line), and the analytical solution for a linear Dμμ (solid line).

32

1

0.8

0.6

0.4 0

0.2

0.4

0.6

0.8

1

Fig. 10. The velocity correlation function normalized to its initial value 3vz (t)vz (0)/v 2 versus the dimensionless time vt/λ . Shown are the functions for the isotropic Dμμ (dotted line), the constant Dμμ (dashed line), and the linear Dμμ (solid line).

33

1 0.8 0.6 0.4 0.2 0 0

2

4

6

8

10

Fig. 11. The running diffusion coefficient versus dimensionless time τ . Shown are the numerical solution for the constant Dμμ (dashed line), the corresponding analytical result as given by Eq. (56) (dotted line), the numerical solution for the linear Dμμ (solid line), and the corresponding analytical result as given by Eq. (107) (dash– dotted line). The units used here are z˜ = zD/v and τ = Dt.

34