Transport of energetic particles in turbulent space plasmas: pitch-angle scattering, telegraph, and diffusion equations

Abstract

Introduction: In this article, we revisit the pitch-angle scattering equation describing the propagation of energetic particles through magnetized plasma. In this case, solar energetic particles and cosmic rays interact with magnetohydrodynamic turbulence and experience stochastic changes in the pitch-angle. Since this happens over an extended period of time, a pitch-angle isotropization process occurs, leading to parallel spatial diffusion. This process is described well by the pitch-angle scattering equation. However, the latter equation is difficult to solve analytically even when considering special cases for the scattering coefficient.

Methods: In the past, a so-called subspace approximation was proposed, which has important applications in the theory of perpendicular diffusion. Alternatively, an approach based on the telegraph equation (also known as telegrapher’s equation) has been developed. We show that two-dimensional subspace approximation and the description based on the telegraph equation are equivalent. However, it is also shown that the obtained distribution functions contain artifacts and inaccuracies that cannot be found in the numerical solution to the problem. Therefore, an N-dimensional subspace approximation is proposed corresponding to a semi-analytical/semi-numerical approach. This is a useful alternative compared to standard numerical solvers.

Results and Discussion: Depending on the application, the N-dimensional subspace approximation can be orders of magnitude faster. Furthermore, the method can easily be modified so that it can be used for any pitch-angle scattering equation.

1 Introduction

The motion of energetic particles such as cosmic rays through plasma is a complicated stochastic process. It is described via transport equations containing different diffusion parameters. The simplest form of a transport equation which is used in this field is the pitch-angle scattering equation ( Shalchi, 2009; Zank, 2014)where we have used time t, particle position along the mean magnetic field z, pitch-angle cosine μ, particle speed v, and pitch-angle scattering coefficient D_μμ. The analytical form of the latter parameter is difficult to determine since it contains information about the interaction between magnetohydrodynamic turbulence and energetic and electrically charged particles. Very originally, a quasi-linear approach was developed to determine the coefficient D_μμ (Jokipii 1966). However, this approach is inaccurate, and it fails to describe correctly the scattering of particles at 90° corresponding to μ = 0 (Shalchi 2009). Therefore, the so-called second-order quasi-linear theory (SOQLT) was developed by Shalchi (2005), which provides non-vanishing scattering at μ = 0, resolving the 90°-problem. This theory was further explored analytically in Shalchi et al. (2009), and the so-called isotropic formwas derived in the limit of a stronger turbulent magnetic field. In Eq. 2, the parameter D does not depend on μ, but it is a complicated function of turbulence and particle properties (Shalchi et al., 2009).

In addition to the question of what the correct analytical form of D_μμ is, one desires to find solutions to Eq. 1. However, so far, no exact solution to the pitch-angle scattering equation has been found, and one has to rely on either a numerical approach or approximations. However, one can show that in the late-time limit, the pitch-angle-averaged distribution functionsatisfies a diffusion or heat transfer equation of the formwhere the parallel spatial diffusion coefficient is related to the pitch-angle scattering coefficient via (Earl, 1974)The heat transfer equation shown above can easily be solved. For sharp initial conditions, for instance, the solution is simply a normalized Gaussian distributioncentered at z = 0 and having the second moment ⟨z²⟩ = 2κ_‖t. One can also write down the more general solutionwhich depends on the initial distribution M(z′, t = 0) and has a Gaussian integral kernel.

More recently (Tautz and Lerche, 2016 and references therein), it was argued that the diffusive solution does not always provide a good approximation, and one should instead use a telegraph equation of the formwhere we have used the telegraph time scaleτ. It should be noted that using the telegraph equation instead of the diffusion equation was, in particular, suggested in the context of adiabatic focusing (Litvinenko and Schlickeiser, 2013; Effenberger and Litvinenko, 2014), but this effect is omitted in this paper.

Independently, a two-dimensional subspace approximation to the solution of Eq. 1 has been developed (see Shalchi et al. (2011) for the original description of this approach and Shalchi (2020) for a review). Although this approach provides only an approximation to the solution of Eq. 1 for the isotropic case, it provides a pitch-angle-dependent solution. The two-dimensional subspace approximation was successfully applied in the theory of perpendicular transport and contributed significantly to the development of advanced particle transport theories (Shalchi, 2020; Shalchi, 2021).

2 The two-dimensional subspace approximation

The two-dimensional subspace approximation was originally developed by Shalchi et al. (2011) to solve pitch-angle scattering Equation 1. We summarize the corresponding results, rewrite previously found solutions, and discuss the relation to the telegraph equation as follows. The following three subsections were mostly taken from Shalchi (2020) but have been modified significantly.

2.1 The isotropic scattering coefficient

For the isotropic scattering coefficient, as given by Eq. 2, the parallel spatial diffusion coefficient is obtained via Eq. 5. Alternatively, we can compute the parallel mean free path that is defined via λ_‖ = 3κ_‖/v. For the isotropic case, those parameters are given byEq. 1 corresponds to a partial differential equation with the variables t, z, and μ. As a first step toward a solution, we use the Fourier transformso that the pitch-angle scattering equation becomesThe inverse Fourier transform is then given by the following equation:For the isotropic scattering coefficient, Eq. 11 is simplified toTo continue, we expand the solution of Eq. 13 in a series of Legendre polynomialswhere the coefficients C_n are functions of time, though they also depend on k. This dependence is not explicitly written down during the following investigations. Using Eq. 14 in the differential Equation 13 yieldswhere denotes the time derivative of the coefficient C_n. In order to further valuate Eq. 15, we use the following two relations for Legendre polynomials (Abramowitz and Stegun, 1974)andWith those two relations, Eq. 15 can be written as follows:To continue, we multiply this equation by the Legendre polynomial P_m, integrate over μ, and use the orthogonality relation of Legendre polynomials (Abramowitz and Stegun, 1974)After performing those steps, we derive the recurrence relationAlternatively, one can use the coefficient Q_m defined viaWith this, the recurrence relation can be written as followsFor the case of no scattering D = 0, we can compare this with the relation (see Equation of Abramowitz and Stegun (1974))where we have used spherical Bessel functions. Thus, we find Q_m = j_m(vkt) for the scatter-free case andUsing this in Eq. 14 yieldswhere we have used Equation 92 from Shalchi et al. (2011). This is also known as plane wave expansion widely used in quantum mechanics. It should be noted that Eq. 25 corresponds to the unperturbed or scatter-free solution. It can be easily obtained directly from Eq. 11 for the case D = 0.

We use Eq. 20 which corresponds to an infinite set of coupled ordinary differential equations. For m = 0, for instance, we findand for m = 1, we obtainIt is problematic here that the coefficients C₀ and C₁ are coupled to C₂. Therefore, it is not possible to derive an exact solution for the coefficients C_n.

2.2 The two-dimensional approximation

Since an exact solution to Eq. 20 seems impossible to be found, one needs to rely on approximations. In the following, we discuss the two-dimensional (2D) subspace approximation originally developed by Shalchi et al. (2011), meaning we setso that only the coefficients C₀ and C₁ are used. In Lasuik and Shalchi (2019), one can find the solution obtained by using a three-dimensional subspace approximation. It is shown that the three-dimensional (3D) solution is too complicated for most applications.

Within the two-dimensional subspace approximation, the expansion (14) is reduced toIn this case, Eqs 26, 27 can be combined to eliminate C₁. Since we set C₂ = 0, we found the second-order differential equationUsing the ansatzleads to the quadratic equationAlternatively, for C₂ = 0, Eqs 26, 27 can be written as the matrix equationAfter using Eq. 31 for both functions C₀(t) and C₁(t), the problem of finding the two ω is expressed asleading to the same quadratic equation as given by Eq. 32. The latter equation can easily be solved by the following equation:We conclude that the eigenvalues can be complex depending on the wave number k. With this, the coefficient C₀ can be written as the linear combinationwith the two unknown coefficients b_±. It follows from Eq. 26 thatThe coefficients b_± will be determined below. Before we perform this task, we write down the solution for F_k(μ, t). We need to combine Eq. 29 with Eqs 36, 37 to deriveIn order to find the coefficients b_±, we can use the initial conditionmeaning that the particle has its initial position at z = 0 and the initial pitch-angle cosine μ₀. Using this in the inverse Fourier transform given by Eq. 12, yields after some straightforward algebraThe latter initial condition used in expansion (Eq. 14) allows us to writeIn order to determine the coefficients C_n(t = 0), we multiply this by P_m and integrate over μ to getTo perform this task, we have used again the orthogonality relation (Eq. 19). For m = 0 and m = 1, this yields¹andTo determine the coefficients b_±, we write down Eqs 36, 37 for t = 0 and use Eqs 43, 44 to deduceThis system of two equations is solved by the following equationUsing this result and Eq. 35 in Eq. 38 provides the two-dimensional subspace approximation to the solution F_k(μ, t). In Section 2.4, we provide a more detailed discussion of this solution.

Our solution is based on the expansion given by Eq. (29). One can easily demonstrate using Eq. 3 andtogether with the orthogonality relation (Eq. 19), that the function C₀(t) corresponds to the Fourier transform of the pitch-angle-averaged distribution function M(z, t), and C₁(t) corresponds to the Fourier transform of the current density or diffusion flux J(z, t). Those two quantities are related to each other via the one-dimensional continuity equationwhich is obtained by averaging Eq. 1 over all μ and using Eqs 3, 47. The exact relations to the coefficients areand

It should be noted that the latter relation is obtained by combining Eq. 47 with Eq. 10 and the expansion given by Eq. 14. After combining these three relations and using the orthogonality relation (Eq. 19), one can obtain Eq. 50. As demonstrated, the coefficients C₀(t) and C₁(t) are directly linked to physical quantities. In particular, the coefficient C₀(t) is very important because it is simply the Fourier transform of the pitch-angle-averaged distribution function M(z, t).

2.3 Further physical quantities

An important quantity in particle transport theory is the characteristic function ⟨e^±ikz⟩. We define the ensemble average viaIt should be noted that in some cases, one could aim for a result that depends on μ₀. Then, the corresponding average is omitted.

To determine the characteristic function, we average over all quantities, and thus, we haveReplacing f(z, μ, t) therein by using Eq. 12 leads toWe now replace F_k(μ, t) by using Eq. 14 and use P₀(μ) = 1 to getDue to the orthogonality of Legendre polynomials (Eq. 19), this is reduced to

To solve the remaining integral, we use Eq. 36 to writeIn order to replace b_±, we use Eq. 46. The integrals over the terms containing μ₀ vanish, and we finally findIt should be noted that the parameters ω_± are given by Eq. 35. For the case that the ω_± are real, the characteristic function given by Eq. 57 is real as well. For the case that the ω_± are complex, it follows from Eq. 35 that . Therefore, the characteristic function is always real, and we have ⟨e^+ikz⟩ = ⟨e^−ikz⟩.

Based on Eq. 35, it can be shown that Eq. 57 contains two asymptotic limits, namely, (see Shalchi (2020) for more details)For small wave numbers, we find the characteristic function of diffusion Equation 4. The result obtained for large wave numbers can be understood as a damped unperturbed orbit.

By comparing Eq. 53 with Eq. 10 and using Eq. 3, we can relate the characteristic function to the μ- and μ₀-averaged functions M(z, t). This relation is given by the following equation:Furthermore, we can compare this with Eq. 49 to findAs an example, we consider the limit D → ∞ so that we can use the first line of Eq. 58 in Eq. 59. We can easily derivecorresponding to a Gaussian solution. The result obtained here is the diffusive solution that one would expect in this case (see Eq. 6 in this paper).

Other physical quantities can be derived by using the subspace approximation, alternative approximations, or even exact calculations (Shalchi, 2006; Shalchi, 2011).

2.4 Rewriting the solution

Eq. 29 corresponds to an integral representation of the solution of the Fourier-transformed pitch-angle scattering equation. This result is based on the 2D subspace approximation. Using therein Eqs 36 and 37, as well as (Eq. 46) yieldswhere the functions ω_± are given by Eq. 35.

The μ- and μ₀-averaged solution is thenin Fourier space. To find the solution in the configuration space, we use Eq. 10 to deriveAlternatively, this result can also be obtained by combining Eqs 57, 59.

It is convenient to define the parameterand it follows from Eq. 35 thatFrom this, we can easily deduceTherewith, the solution in the configuration space is given as the following Fourier transformwhere we have used the integrand, which is an even function of k. With the help of hyperbolic functions, this can be written in a more compact formIt should be noted that the quantity ξ, given by Eq. 65, can either be real or imaginary depending on the value of k. Eq. 69 provides an integral representation of the μ- and μ₀-averaged distribution function based on the 2D subspace approximation. Alternative forms are presented in Supplementary Appendix S1 of this paper. In Supplementary Appendix S2, we provide an approximative solution of the remaining integral.

2.5 Relation to the Telegrapher’s equation

We have derived an ordinary differential equation for the function C₀(t) previously (Eq. 30), which can be written as followsAs shown via Eq. 49, C₀(t) is the Fourier transform of the distribution function M(z, t). Thus, working in the configuration space instead of the Fourier space allows us to write Eq. 70 asThe latter equation has the same form as Eq. 8, and, thus, it corresponds to a telegraph equation. A quick alternative derivation of the latter equation can be found in Supplementary Appendix S3. It should be noted that the coefficient C₀(t) used here depends also on the initial pitch-angle cosine μ₀. If one averages over the latter quantity, the two-dimensional subspace approximation provides Eq. 69. In Supplementary Appendix S4, we demonstrate that the latter form indeed solves Eq. 71. Using therein Eq. 9 and the scattering time τ = 1/(2D) yields the telegraph equation, as given by Eq. 8. Therefore, we have shown the complete equivalence of the two-dimensional subspace approximation and the telegraph equation. The solution given by Eq. 69 is correctly normalized. In order to demonstrate this, we considerTherein, we use (Zwillinger, 2012)to write this asFrom Eq. 65, it follows that ξ_k=0 = D, and, thus, we findAs already pointed out in Tautz and Lerche (2016), one can use the transformationand Eq. 8 becomesThis corresponds to the Klein–Gordon equation but with imaginary mass. After comparing Eqs 69, 76 with each other, we can easily read off the function .

We have focused on the function C₀(t) previously. We can also derive an ordinary differential equation for C₁(t). By combining Eqs 26, 27, we derivewhere we have set C₂ = 0 corresponding to the 2D subspace approximation. Eq. 78 is the same equation as we have derived above for C₀. The function C₁(t) corresponds to the Fourier-transformed current density, as shown by Eq. 50. Therefore, the telegraph and Klein–Gordon equations can also be derived for the current density function. In order to obtain the current density, as a further solution to the telegraph equation, we can combine Eq. 69 with the continuity Equation 48. We can easily derivewhere ξ is given by Eq. 65. Of course, integrating the obtained over all z yields 0, meaning that the found solution to the telegraph equation is not normalized to 1.

3 The N-dimensional subspace approximation

Previously, we have used the expansion in the Legendre polynomials (see Eq. 14 of this paper). The functions C_n(t) therein are given by the recurrence relation (Eq. 20). This relation is still exact and can be written as the following matrix equationwith the tridiagonal matrix A having the componentsThe vector in Eq. 80 contains the functions C_n(t) needed in the expansion given by Eq. 14. The formal solution of Eq. 80 can be written as follows:where we have used the matrix exponential. The initial conditions C_n(t = 0) are given by Eq. 42. Eq. 82 can be easily evaluated with software such as MATLAB. However, it is required to work with a finite matrix A. This corresponds to the subspace approximation outlined above. Let us assume that we work with an N × N-matrix. This then corresponds to an N-dimensional subspace approximation. The method described here corresponds to a semi-numerical/semi-analytical approach that solves the pitch-angle scattering equation, but this method can be faster if one needs the solution only for a specific time t. Standard numerical solvers (see Supplementary Appendix S5 of this paper) require a high time resolution to be accurate. Therefore, one typically needs to work with roughly thousand time-steps so that the solution converges to the true solution of the differential equation. The N-dimensional subspace approximation can be applied to a single time value. As shown via Figures 2–9, an accurate solution is obtained for N = 10.

For certain applications, one could be interested in the μ- and μ₀-averaged solution only. Analytical solutions of diffusion and telegraph equations are incomplete and inaccurate depending on the considered application. For the case of pitch-angle-averaged solutions, the N-dimensional subspace approximation is particularly powerful, as outlined below. First, we define the matrix exponential used already above viaThen, Eq. 82 can be written as followsor in component notation,At the initial time, the components of the vector are given by Eq. 42. If those coefficients are averaged over μ₀, we can easily derivemeaning that all coefficients are 0, except C₀(t = 0). Therefore, we can write the time-dependent coefficients asFurthermore, the μ-dependent solution is given by Eq. 14. After μ-averaging of the latter expansion, we findwhere M_k(t) is the Fourier-transformed distribution function as observed by Eq. 49. It should be noted that the function C₀(t) discussed here is also μ₀-averaged. Furthermore, the characteristic function is easily obtained viameaning that the 00-component is simply the characteristic function. Thus, it follows from Eq. 59 thatwhich is an integral and matrix exponential representation of the μ- and μ₀-averaged distribution function. Therefore, in order to obtain the distribution function M(z, t) for given z and t, we need to numerically solve the k-integral in Eq. 90. For each value of k, we set up the matrix A defined via Eq. 81, numerically compute the matrix exponential E, and use the component E₀₀ in the numerically evaluated k-integral. The distribution functions shown in Figures 6, 7, based on the 10D subspace approximation, for instance, can be computed with a regular computer within a few seconds. Figure 1 shows a comparison in speed between different numerical methods. This comparison includes the N-dimensional subspace approximation described above and the pure numerical approach described in Supplementary Appendix S5 of this paper, which corresponds to an implicit Euler method.

FIGURE 1

The μ- and μ₀-dependent Fourier-transformed solution is given by the following equation:where we have combined Eqs 14, 42, and 85. The Fourier transform can be performed using Eq. 10 and solving the k-integral numerically. Of course, obtaining and plotting the pitch-angle-dependent result is more time-consuming when compared to the pitch-angle-averaged solution.

In certain analytical theories developed for describing the perpendicular transport of energetic particles, one needs to know the function (Shalchi, 2010; Shalchi, 2017; Shalchi, 2020; Shalchi, 2021)that is somewhat similar but not identical compared to the characteristic function discussed above. In order to express Γ_k(t) as before, we perform the same mathematical steps. The pitch-angle-dependent solution is given by Eq. 14. In order to obtain Γ_k(t), we needTo evaluate this further, we use Eqs 12, 14. After using those two relations, we deriveFor the μ-integral, we can use the orthogonality relation (Eq. 19) to findUsing the above relation allows us to perform the following steps:where we have used Eqs 19, 42, and 85. Therefore, the derived function corresponds to the matrix element E₁₁ which can be computed quickly based on the N-dimensional subspace approximation. Figure 12 shows some example plots for the quantity .

4 Comparison of results

We have solved the pitch-angle scattering equation numerically using an implicit Euler method (Supplementary Appendix S5) and the N-dimensional subspace approximation outlined in the previous section. We have considered two cases, namely, N = 2 (corresponding to the pure analytical case discussed above) and N = 10 (which provides an accurate result). In most cases, we have only considered the μ-averaged solution to reduce the number of plots. Some results are also averaged over the initial pitch-angle cosine μ₀.

Figures 2–5 show the characteristic function ⟨e^ikz⟩≡⟨e^−ikz⟩, which corresponds to the Fourier transform of the distribution function . In Figures 2, 3, the characteristic function is plotted versus the dimensionless wave number for different values of the dimensionless time . We have also shown the solution of the diffusion equation as given by the first line of Eq. 58. We can easily see that all solutions agree with each other at later times. This is not the case for early times where the 2D subspace approximation and the diffusive solution differ significantly from the numerical solution. The 10D subspace approximation agrees very well with the numerical solution in all considered cases.

Figures 4, 5 show the characteristic function versus time for different values of . We can easily see agreement for smaller values of but disagreement for larger values. However, the 10D subspace approximation and the numerical solution agree very well with each other. It has to be emphasized that the 10D subspace approximation solution, which can be seen as a semi-analytical/semi-numerical technique, is several orders of magnitude faster than standard numerical solvers.

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

Figures 6–9 show the distribution function . Figures 6, 7 show this function versus the dimensionless position for different times. For late times, all considered results agree with each other, as expected. The corresponding distributions are well-described by the Gaussian given by Eq. 6. For early times, however, diffusive and 2D subspace results do not agree well with the pure numerical solution. The 2D subspace solution contains spikes that are a consequence of the Dirac delta (Supplementary Appendix S2). The diffusive solution is non-zero everywhere in space. Numerical and 10D subspace solutions correctly describe that the distribution function is exactly 0 for due to the finite propagation speed of the particles. The latter effect can be observed much better by plotting the distribution function versus time for different values of . This is done via Figures 8, 9.

FIGURE 6

FIGURE 7

FIGURE 8

FIGURE 9

Figure 10 shows the comparison of the time evolution of based on diffusion equation and the 10D subspace approximation. We can clearly see the similarity for later times. For early times, on the other hand, we observe significant differences. In particular, the 10D solution provides as needed. Alternatively, we have plotted versus time for different positions (Figure 11) where the aforementioned effect can be observed more clearly.

FIGURE 10

FIGURE 11

Last but not the least, we have computed the function defined via Eq. 92. The latter function enters certain analytical theories for perpendicular diffusion. According to Figure 12, the 2D subspace approximation works overall well for computing this quantity. This explains why analytical theories for perpendicular diffusion, in which the 2D subspace approximation was used, agree well with performed test-particle simulations (Shalchi, 2020; Shalchi, 2021).

FIGURE 12

5 Summary and conclusion

In this paper, we have focused on the most basic transport equation, namely, the pitch-angle scattering equation, as given by Eq. 1. Analytical and numerical investigations of pitch-angle-dependent transport equations have been the subject of several papers published during recent years. In addition to studies of the basic pitch-angle scattering equation (Shalchi et al., 2011; Tautz and Lerche, 2016; Lasuik and Shalchi, 2017; Lasuik and Shalchi, 2019), authors have explored the impact of so-called focusing, an effect which is related to a non-constant mean magnetic field (Danos et al., 2013; Litvinenko and Schlickeiser, 2013; Effenberger and Litvinenko, 2014; Lasuik et al., 2017; Wang and Qin, 2020; Wang and Qin, 2021; Wang and Qin, 2023). Even more complicated cases, including perpendicular particle transport, have been investigated by Wang and Qin (2024).

In this article, we have reviewed the two-dimensional subspace approximation originally developed by Shalchi et al. (2011) and discussed the provided solutions in configuration and Fourier spaces. We have also demonstrated that the two-dimensional subspace approximation is equivalent to using a telegraph equation. Normalized solutions in configuration and Fourier spaces are also discussed. However, we also argue that such solutions do not often provide appropriate results even if the pitch-angle average is considered. Although it was often argued that the telegraph equation is more complete than the usual diffusion approach (Tautz and Lerche, 2016), the solution discussed here contains artifacts that are not realistic. In particular, we observe spikes at z = ±vt (Figures 6, 7 as well as Supplementary Appendix S2).

Therefore, it is important to solve the pitch-angle scattering equation numerically. However, standard approaches such as implicit Euler or Crank–Nicolson solvers are time-consuming to use. In this paper, we have, thus, developed an N-dimensional subspace approach. This method can be seen as a semi-analytical/semi-numerical method. It has the advantage of being is several orders of magnitude faster than standard solvers (see Figure 1 of this paper). This is in particular the case if one is only interested in pitch-angle-averaged solutions at a given time. Standard solvers require a high time resolution in order to provide an accurate result. The N-dimensional subspace technique can be applied for a single time value if this is everything what is needed. It should also be emphasized that the N-dimensional subspace method can be easily parallelized since for a given k and t, one can compute the matrix exponentials independently of other values. This is also valid if one is looking for a μ- and μ₀-dependent result.

In this paper, we have computed several quantities such as distribution and characteristic functions as well as the function which is defined via Eq. 92 of this paper. We have compared numerical solutions obtained by using a standard solver with results obtained by using the N-dimensional subspace approximation for different values of N and the diffusive solution. The main difference is that pure numerical and 10D solutions provide , meaning that the particles have a finite propagation speed. All results are visualized in Figures 2–12. One can clearly see that for N = 10, we obtain an accurate result that agrees well with the pure numerical solution of the pitch-angle scattering equation.

It has to be noted that the N-dimensional subspace method presented in this paper was specifically developed for the basic pitch-angle scattering equation and an isotropic scattering coefficient. However, this approach can be easily modified so that it can be used for more general transport equations including focused transport equations and other forms of the pitch-angle scattering coefficient.

References

• 1

  AbramowitzM.StegunI. A. (1974). Handbook of mathematical functions. New York: Dover Publications.

  • Google Scholar

• 2

  CrankJ.NicolsonP. (1947). A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Camb. Phil. Soc.43, 50–67. 10.1017/s0305004100023197

  • CrossRef
  • Google Scholar

• 3

  DanosR. J.FiegeJ. D.ShalchiA. (2013). Numerical analysis of the fokker-planck equation with adiabatic focusing: isotropic pitch-angle scattering. Astrophysical J.772, 35. 10.1088/0004-637x/772/1/35

  • CrossRef
  • Google Scholar

• 4

  EarlJ. A. (1974). The diffusive idealization of charged-particle transport in random magnetic fields. Astrophys. J.193, 231. 10.1086/153152

  • CrossRef
  • Google Scholar

• 5

  EffenbergerF.LitvinenkoY. E. (2014). The diffusion approximation versus the telegraph equation for modeling solar energetic particle transport with adiabatic focusing, I. Isotropic pitch-angle scattering. Astrophysical J.783, 15. 10.1088/0004-637x/783/1/15

  • CrossRef
  • Google Scholar

• 6

  GradshteynI. S.RyzhikI. M. (2000). Table of integrals, series, and products. New York: Academic Press.

  • Google Scholar

• 7

  JokipiiJ. R. (1966). Cosmic-Ray propagation. I. Charged particles in a random magnetic field. Astrophys. J.146, 480. 10.1086/148912

  • CrossRef
  • Google Scholar

• 8

  LasuikJ.FiegeJ. D.ShalchiA. (2017). Numerical analysis of the Fokker-Planck equation with adiabatic focusing: realistic pitch-angle scattering. Adv. Space Res.59, 722–735. 10.1016/j.asr.2016.10.027

  • CrossRef
  • Google Scholar

• 9

  LasuikJ.ShalchiA. (2017). Solutions of the cosmic ray velocity diffusion equation. Adv. Space Res.60, 1532–1546. 10.1016/j.asr.2017.06.035

  • CrossRef
  • Google Scholar

• 10

  LasuikJ.ShalchiA. (2019). Subspace approximations to the cosmic ray Fokker-Planck equation. Mon. Notices R. Astronomical Soc.485, 1635–1650. 10.1093/mnras/stz474

  • CrossRef
  • Google Scholar

• 11

  LitvinenkoY. E.SchlickeiserR. (2013). The telegraph equation for cosmic-ray transport with weak adiabatic focusing. A A554, A59. 10.1051/0004-6361/201321327

  • CrossRef
  • Google Scholar

• 12

  PressW. H.TeukolskyS. A.VetterlingW. T.FlanneryB. P. (2007). Numerical Recipes. Cambridge University Press.

  • Google Scholar

• 13

  ShalchiA. (2005). Second-order quasilinear theory of cosmic ray transport. Phys. Plasmas12, 052905. 10.1063/1.1895805

  • CrossRef
  • Google Scholar

• 14

  ShalchiA. (2006). Analytical investigation of the two-dimensional cosmic ray Fokker-Planck equation. Astronomy Astrophysics448, 809–816. 10.1051/0004-6361:20053664

  • CrossRef
  • Google Scholar

• 15

  ShalchiA. (2009). Nonlinear cosmic ray diffusion theories. Astrophysics Space Sci. Libr.362. Springer, Berlin. 10.1007/978-3-642-00309-7

  • CrossRef
  • Google Scholar

• 16

  ShalchiA. (2010). A unified particle diffusion theory for cross-field scattering: subdiffusion, recovery of diffusion, and diffusion in 3D turbulence. Astrophys. J.720, L127–L130. 10.1088/2041-8205/720/2/l127

  • CrossRef
  • Google Scholar

• 17

  ShalchiA. (2011). Velocity correlation functions of charged particles derived from the fokker-planck equation. Adv. Space Res.47, 1147–1164. 10.1016/j.asr.2010.12.002

  • CrossRef
  • Google Scholar

• 18

  ShalchiA. (2017). Time-dependent perpendicular transport of energetic particles in magnetic turbulence with transverse complexity. Phys. Plasmas24, 050702. 10.1063/1.4982805

  • CrossRef
  • Google Scholar

• 19

  ShalchiA. (2020). Perpendicular transport of energetic particles in magnetic turbulence. Space Sci. Rev.216, 23. 10.1007/s11214-020-0644-4

  • CrossRef
  • Google Scholar

• 20

  ShalchiA. (2021). Perpendicular diffusion of energetic particles: a complete analytical theory. Astrophys. J.923, 209. 10.3847/1538-4357/ac2363

  • CrossRef
  • Google Scholar

• 21

  ShalchiA.ŠkodaT.TautzR. C.SchlickeiserR. (2009). Analytical description of nonlinear cosmic ray scattering: isotropic and quasilinear regimes of pitch-angle diffusion. Astron. Astrophys.507, 589–597. 10.1051/0004-6361/200912755

  • CrossRef
  • Google Scholar

• 22

  ShalchiA.TautzR. C.RempelT. J. (2011). Test-particle transport: higher-order correlations and time-dependent diffusion. Plasma Phys. control. Fusion53, 105016. 10.1088/0741-3335/53/10/105016

  • CrossRef
  • Google Scholar

• 23

  TautzR. C.LercheI. (2016). Application of the three-dimensional telegraph equation to cosmic-ray transport. Res. Astronomy Astrophysics16, 162. 10.1088/1674-4527/16/10/162

  • CrossRef
  • Google Scholar

• 24

  WangJ. F.QinG. (2020). The invariance of the diffusion coefficient with iterative operations of the charged particle transport equation. Astrophysical J.899, 39. 10.3847/1538-4357/aba3c8

  • CrossRef
  • Google Scholar

• 25

  WangJ. F.QinG. (2021). Study of momentum diffusion with the effect of adiabatic focusing. Astrophysical J. Suppl. Ser.257, 44. 10.3847/1538-4365/ac1bb3

  • CrossRef
  • Google Scholar

• 26

  WangJ. F.QinG. (2023). Relationship of transport coefficients with statistical quantities of charged particles. Astrophysical J.954, 213. 10.3847/1538-4357/ace9d3

  • CrossRef
  • Google Scholar

• 27

  WangJ. F.QinG. (2024). The effect of solar wind on charged particles’ diffusion coefficients. Astrophysical J.961, 6. 10.3847/1538-4357/ad09b7

  • CrossRef
  • Google Scholar

• 28

  WatsonG. N. (2011). A treatise on the theory of Bessel functions. Cambridge: Cambridge Mathematical Library.

  • Google Scholar

• 29

  ZankG. P. (2014). Transport processes in space Physics and astrophysics. Lect. Notes Phys.877. Springer, New York. 10.1007/978-1-4614-8480-6

  • CrossRef
  • Google Scholar

• 30

  ZwillingerD. (2012). Standard mathematical tables and formulae. Boca Raton: CRC Press.

  • Google Scholar