Transport in Astrophysics: I. Diffusion of Solar and Galactic Cosmic Rays

Lorenzo Zaninetti

Physics Department, Turin, Italy.

DOI: 10.4236/ijaa.2022.121003                

Abstract

Some solutions for the diffusion phenomenon as a function of time and space are reviewed. Two new solutions of the homogeneous diffusion equation in 1D and 2D are derived in the presence of an existing fixed number of particles. The initial conditions which allow deriving a power law behavior for the energy of the cosmic rays (CR) are derived. The superposition of transient diffusive phenomena on an existing power law distribution for the energy of CR allows simulating the knee, the second knee, and the ankle.

Keywords

Cosmic Rays, Particle Diffusion, Random Walks

1. Introduction

We will briefly review the meaning of typical features of Cosmic Rays (CR). After their discovery in 1912 by Hess [1], this term started to appear in [2] [3] [4] [5]. The term “solar modulation” is connected with the solar activity and its 11-year cycle; this term was introduced by [6] followed by [7] [8] and an example is visible in Figure 6 in [9]. The solar modulation of low energy CR, as an example 0.32 GeV, can decrease the intensity by a factor of ≈20, see Figure 5 in [10]. The term “knee” in the CR energy distribution was introduced for the first time in 1961 by [11] and subsequently widely used, see as an example [12] [13] [14] [15]. The term “ankle” refers to a change in the energy distribution of CR at ≈10⁷ TeV; this term was introduced in 1996 by [16] followed by [17] [18]. The term “galactic cosmic rays” was introduced by Compton [19] in order to predict asymmetries in the intensities of CR due to the motion of our sun with a velocity of 300 km/sec, which now is taken to be 230 km/s; this trend continued with [20], who introduced the scattering due to the stars and with [21] [22] where the radio-emission of our galaxy radiation was supposed to be due to the gyration of CR around magnetic fields. The term “extra-galactic cosmic rays” was introduced by Burbidge [23] where CR with energies between 10¹⁸ eV and 10²⁰ eV are supposed to be accelerated in the clusters of galaxies; [24] estimated the rate of CR production by extra-galactic sources as well as by active galaxies, [25] suggested that CR were accelerated by galactic and extra-galactic supernovae. The first models for the acceleration of CR were assumed to be the encounters of relativistic particles with astrophysical clouds [26], followed by the presence of plasma oscillations [27], the interaction with hydro-magnetic waves [28] and the resonance with magneto-hydrodynamic waves [29]. The term “diffusion of cosmic rays” was introduced by [30], where bursts of CR occurring randomly in time instead of continuously were adopted. The research on diffusion then continued; we select some approaches, among others: the statistical homogeneity and isotropy of the magnetic fluctuations were analysed in [31], the anisotropic diffusion of CR in the interplanetary medium due to the irregular spiral interplanetary magnetic field [32], the ordinary diffusion tensor [33], the presence of a magnetic barrier in interplanetary space [34], the existence of magnetic fields from meteor streams [35], an anisotropic-diffusion approximation [36], 3D random walk of the interstellar magnetic field lines [37], an analytic expression for the power spectrum, $P\left(k\right)$, [38], the modulation of CR by an interplanetary shock wave [39], time profiles in the intensity of solar CR at various distances from the source as a function of the ratio of the mean free path to the focusing length of the interplanetary field [40], solution of the diffusion equations adopting realistic models of the galactic field and using diffusion coefficients appropriate for strong turbulence [41], the computation of the perpendicular diffusion coefficients and mean free paths of particles for an anisotropic Alfvénic turbulence spectrum, [42], asymmetric diffusion in the presence of high-amplitude magneto-hydrodynamic turbulence [43], a two-component model for the evolution of fluctuations of solar wind plasma [44] and anomalous transport phenomena associated with galactic CR propagating through interstellar space [45]. The present paper reviews the existing situation of the solutions of the diffusion equation and derives two new solutions in 1D and 2D, see Section 2. The astrophysical applications to CR allow building an energy spectrum similar to the observed one, and simulate the knee, the second knee and the ankle, see Section 3.

2. Transient Diffusion

This section reviews the definition of the diffusion coefficient, Fick’s second law for diffusion, the impulsive 1D, 2D and 3D solutions, introduces two new solutions (2D and 3D) for the impulsive case in the presence of an existing fixed profile and reviews the diffusion internal to a ball when the density on the boundary is constant over time.

2.1. The Diffusion Coefficient

The dependence for the mean square displacement, $\overline{R^2\left(t\right)}$, according to Equation (8.38) in [46] is

$\overline{R^2\left(t\right)}=2dDt\left(t\to \infty \right)\mathrm{,}$ (1)

where d is the number of spatial dimensions. From Equation (1), the diffusion coefficient is derived in the continuum:

$D=\left(t\to \infty \right)\frac{{\bar{R}}^2}{2dt}\mathrm{.}$ (2)

Using discrete time steps, the average square radius after N steps, Equation (12.5) in [46], is

$\langle R^2\left(N\right)\rangle \sim 2dDN\mathrm{,}$ (3)

from which the diffusion coefficient is derived:

$D=\frac{\langle R^2\left(N\right)\rangle }{2dN}.$ (4)

If $\langle R^2\left(N\right)\rangle \sim \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}N$, the diffusion coefficient is

$D=\frac{1}{2d}\lambda v_{tr}\mathrm{,}$ (5)

when the step length of the walker or mean free path between successive collisions is $\lambda$ and the transport velocity is $v_{tr}$.

2.2. Fick’s Second Law

Fick’s second law in 3D states that a change in concentration, N, in any part of the system is due to an inflow and an outflow of material into and out of that part of the system

$\begin{array}{l}\frac{\partial }{\partial t}N\left(x,y,z,t\right)=D{\nabla }^{2}N\\ =D\left(\frac{{\partial }^2}{\partial x^2}N\left(x,y,z,t\right)+\frac{{\partial }^2}{\partial y^2}N\left(x,y,z,t\right)+\frac{{\partial }^2}{\partial z^2}N\left(x,y,z,t\right)\right),\end{array}$ (6)

where D is the diffusion coefficient, t is the time and ${\nabla }^2$ is the Laplacian operator.

2.3. 1D Case, Impulsive Injection

In 1D, Fick’s second law is

$\frac{\partial }{\partial t}N\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)=D\left(\frac{{\partial }^2}{\partial x^2}N\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)\right)\mathrm{,}$ (7)

and a first solution is of Gaussian type

$N\left(x\mathrm{,}t\right)=\frac{N{\text{e}}^{-\frac{x^2}{4Dt}}}{2\sqrt{\pi Dt}}\mathrm{,}$ (8)

where $N=i\times dt$ is the total number of particles injected in the time dt and i the rate of particles injected at the centre.

At position x the concentration has a maximum at $t=t_{\max }$ where

$t_{\max }=\frac{x^2}{2D}.$ (9)

Figure 1 reports the number of particles for different times.

A second solution is given by

$N\left(x\mathrm{,}t\right)=\frac{N_0\text{erfc}\left(\frac{x}{2\sqrt{Dt}}\right)}{2}\mathrm{,}$ (10)

where $N_0$ is the concentration at $N\left(\mathrm{0,0}\right)$, $\text{erfc}\left(x\right)$ is the complementary error function defined by

$\text{erfc}\left(x\right)=1-\text{erf}\left(x\right)$ (11)

and $\text{erf}\left(x\right)$ is the error function [47]. Figure 2 reports a comparison of the two solutions here analysed.

2.4. 2D Case, Impulsive Injection

In 2D, Fick’s second law in polar coordinates is

$D\left(\frac{\frac{\partial }{\partial r}N\left(r\mathrm{,}t\right)}{r}+\frac{{\partial }^2}{\partial r^2}N\left(r\mathrm{,}t\right)\right)=\frac{\partial }{\partial t}N\left(r\mathrm{,}t\right)\mathrm{.}$ (12)

A solution in the impulsive case is

Figure 1. Number of particles as function of the distance (pc) for different times (years) as in the legend when $N=1$ and $D=\mathrm{1\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$.

Figure 2. Number of particles as function of the distance (pc) for the Gaussian solution when $N=1$, $D=\mathrm{1\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$ (green line), and for the error solution solution when $N_0=1$, $D=\mathrm{1\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$ (red line).

$N\left(r\mathrm{,}t\right)=\frac{N_0{\text{e}}^{-\frac{r^2}{4Dt}}}{4\pi Dt}$ (13)

where $N_0$ is the number of particles injected in the time dt. The concentration as a function of r has a maximum at $t=t_{\max }$ where

$t_{\max }=\frac{r^2}{4D}.$ (14)

2.5. 3D Case, Impulsive Injection

In 3D, Fick’s second law in spherical coordinates is

$D\left(\frac{2\left(\frac{\partial }{\partial r}N\left(r\mathrm{,}t\right)\right)}{r}+\frac{{\partial }^2}{\partial r^2}N\left(r\mathrm{,}t\right)\right)=\frac{\partial }{\partial t}N\left(r\mathrm{,}t\right)\mathrm{,}$ (15)

which has a solution

$N\left(r\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}D\right)=\frac{N_0\sqrt{4}{\text{e}}^{-\frac{r^2}{4Dt}}}{16{\left(\pi Dt\right)}^{\frac{3}{2}}}\mathrm{,}$ (16)

where $N_0$ is the number of particles injected at $t=0$ and $r=0$. Once the radius of the sphere, r is fixed, the maximum of the number of particles is at time

$t_{\max }=\frac{r^2}{6D}.$ (17)

2.6. 1D Case, Existing Profile

We now solve Fick’s second law in 1D as given by Equation (7) over the spatial domain $\left[\mathrm{0,}L\right]$ in the presence of an existing profile (the initial condition) in the number of particles

$N\left(x\mathrm{,0}\right)=N_0{\text{e}}^{-\frac{xs}{L}}$ (18)

where s is an adjustable parameter and $N_0$ the number of particles at $x=0$.

The boundary conditions are assumed to be $\frac{\partial }{\partial x}u\left(\mathrm{0,}t\right)=0$ and $\frac{\partial }{\partial x}u\left(L\mathrm{,}t\right)=0$ and the solution is

$u\left(x\mathrm{,}t\right)=-\frac{N_0\left({\text{e}}^{-s}-1\right)}{s}+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(-\frac{2\cos \left(\frac{n\pi x}{L}\right){\text{e}}^{-\frac{D{\pi }^2{n}^{2}t}{L^2}}s{N}_0\left({\text{e}}^{-s}{\left(-1\right)}^n-1\right)}{{\pi }^2{n}^2+s^2}\right)}$. (19)

Figure 3 reports the number of particles for different times.

Figure 3. Number of particles in the case of an existing profile as function of the distance (pc) for different times (years) when $N_0=1$ for different times as in the legend $D=\mathrm{1\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$.

2.7. 2D Case, Existing Profile

We now solve Fick’s second law in 2D, which is

$\frac{\partial }{\partial t}N\left(x,y,t\right)=D\left(\frac{{\partial }^2}{\partial x^2}N\left(x,y,t\right)+\frac{{\partial }^2}{\partial y^2}N\left(x,y,t\right)\right),$ (20)

assuming the following profile of density at $t=0$

$N\left(x,y,0\right)=\frac{x\left(L-x\right)\left(L-y\right)y}{L^4},$ (21)

where L is the side of the square. The boundary conditions are assumed to be $N\left(0,y,t\right)=0$, $N\left(L,y,t\right)=0$, $N\left(x,0,t\right)=0$ and $N\left(x,L,t\right)=0$ and the solution is

$\begin{array}{l}N\left(x\mathrm{,}y\mathrm{,}t\right)\\ ={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(-\frac{16\sin \left(\frac{n\pi x}{L}\right)\sin \left(\frac{m\pi y}{L}\right){\text{e}}^{-\frac{D{\pi }^{2}t\left(m^2+n^2\right)}{L^2}}\left(-{\left(-1\right)}^{m+n}+{\left(-1\right)}^m+{\left(-1\right)}^n-1\right)}{m^3{\pi }^6{n}^3}\right)}}.\end{array}$ (22)

An example is reported in Figure 4.

2.8. 3D Case, in the Ball

The propagation of the temperature in a ball has been investigated, see

https://www.math.hmc.edu/ajb/PCMI/lecture_schedule.html, and due to the analogy between the heat equation and the diffusion equation, we adopt the same

Figure 4. Number of particles as function of x and y when $t=1.15\times {10}^5\text{yr}$, $D={10}^{-2}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$ and $m=n=10$.

solution. The variable are: the radius of the ball, a, the variable radius, r, the initial number of particles in the ball, $N_b$, the boundary number of particles which is constant with time, $N_0$, the diffusion coefficient, D, and the index of the Fourier series, n. The solution to Equation (6) is

$N\left(r,t\right)=N_b+\frac{\left(2N_0-2N_b\right)a\left({\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{{\left(-1\right)}^{n+1}{\text{e}}^{-\frac{D{n}^2{\pi }^{2}t}{a^2}}\sin \left(\frac{n\pi r}{a}\right)}{n}}\right)}{\pi r}.$ (23)

Figure 5 reports an example of the time necessary to fill the ball.

3. Astrophysical Applications

In the following, the space is expressed in pc and the time in years.

3.1. The Relativistic Gyroradius

The relativistic gyroradius or Larmor radius is

$r_H=\frac{mc\gamma v_{\perp }}{qB}\mathrm{,}$ (24)

where m is the mass of the particle, c is the speed of light,

$\gamma =\frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\mathrm{,}$ (25)

Figure 5. Number of particles as function of r andt as reported in the legend when $D={10}^{-1}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$, $N_0=1000$, $N_b=1$, $a=100\text{pc}$, and 200 iterations.

is the Lorentz factor, v is the velocity of the particle,q is the charge of the particle, B is the magnetic field and $v_{\perp }$ is the velocity perpendicular to the magnetic field, see formula (1.54) in [48] or formula (7.3) in [49]. In the case of CR we express the Larmor radius in pc

$r_L=\frac{1.081\times {10}^{-6}{E}_{\text{GeV}}}{Z{B}_{-6}}\text{pc}=\frac{1.081E_{15}}{Z{B}_{-6}}\mathrm{,}$ (26)

where Z is the atomic number, $E_{\text{GeV}}$ is the energy expressed in GeV, $E_{15}$ is the energy expressed in 10¹⁵ eV, and $B_{-6}$ is the magnetic field expressed in 10⁻⁶ gauss, see [50]. On assuming that the CR diffuse with a mean free path equal to the relativistic gyroradius, the transport velocity is equal to the speed of light and $d=3$, the diffusion coefficient according to Equation (5), is

$D=\frac{0.055134E_{15}}{Z{B}_{-6}}\frac{{\text{pc}}^{\text{2}}}{\text{year}}=\frac{5.5134\times {10}^{-8}{E}_{\text{GeV}}}{Z{B}_{-6}}\frac{{\text{pc}}^{\text{2}}}{\text{year}}.$ (27)

3.2. Cosmic Rays

The observed differential spectrum of CR according to [51] [52] [53] [54] is reported in Figure 6 in the H case ( $I_H$ ).

Flux of H versus energy per nucleus in Gev: experimental data (empty stars) and theoretical power law (full line), see Figure 6.

From this fit it is possible to derive an approximate behaviour for the H intensity in the range [5 GeV - 3.93 × 10⁵ GeV]:

$I_H\left(E\right)\approx 1.24\mp {10}^4{\left(\frac{E}{1\text{GeV}}\right)}^{-2.74}\frac{\text{nucleons}}{{\text{m}}^{\text{2}}\cdot \text{s}\cdot \text{GeV}\cdot \text{sr}}\mathrm{.}$ (28)

Figure 6. Flux of H versus energy per nucleus in Gev: experimental data (empty stars) and theoretical power law (full line).

3.3. The Spectral Index

In the 3D impulsive case, see solution (16), the ratio between the two values of concentration characterized by $D_{\min }$ and $D_{\max }$ is

$\frac{N\left(r,t;D_{\min }\right)}{N\left(r,t;D_{\max }\right)}=\frac{{\text{e}}^{-\frac{r^2\left(D_{\max }-D_{\min }\right)}{4D_{\min }t{D}_{\max }}}{D}_{\max }^{\frac{3}{2}}}{D_{\min }^{\frac{3}{2}}}.$ (29)

This ratio can be parametrized as

$\frac{N\left(r,t;D_{\min }\right)}{N\left(r,t;D_{\max }\right)}=\frac{1}{1000},$ (30)

which means a spectral index in the differential number of CR $\alpha =-3$. The solution of the above equation when $D_{\min }=0.02$, $D_{\max }=0.2$ and $N_0=1$ is

$t=1.0857r^2,$ (31)

and Figure 7 reports this relation. The behaviour of the spectral index as a function of the time can be easily evaluated:

$\alpha ={\log }_{10}\left(\frac{N\left(r,t;D_{\min }\right)}{N\left(r,t;D_{\max }\right)}\right),$ (32)

see Figure 8.

This figure shows that diffusion alone does not produce a stable spectral index of CR.

Figure 7. The locus of time versus distance, r, for which the 3D diffusion from a point gives $I_H\left(E\right)\propto {\left(\frac{E}{1\text{GeV}}\right)}^{-3}$ when $D_{\min }=0.02$, $D_{\max }=0.2$ and $N_0=1$.

Figure 8. The spectral index which characterizes the impulsive 3D diffusion from a point, $I_H\left(E\right)\propto {\left(\frac{E}{1\text{GeV}}\right)}^{\alpha }$, when $r=10\text{pc}$, $D_{\min }=0.02$ and $D_{\max }=0.2$.

3.4. Modulation of CR at Low Energies

The spectrum of CR without the influence of the heliosphere was measured in the summer of 2012 by the Voyager 1 spacecraft, see [55] [56] [57], and is reported in Figure 9 where the range [2.06 × 10⁻³ GeV - 0.58 GeV] is covered. In order to have a continuous relation between the number of existing CR, $N_{exi}$, and the energy, this observational data has been fitted by a polynomial regression [58]:

$N_{exi}=a_1+a_{2}E+a_3{E}^2+a_4{E}^3+a_5{E}^4,$ (33)

where $a_1=3.541$, $a_2=-1.135$, $a_3=-0.278$, $a_4=0.122$ and $a_5=3.992\times {10}^{-2}$, see Figure 10.

In the framework of the 3D impulsive case, see solution (16), we now report in Figure 11 the number of CR at 1 au, $N_{imp}$, as a function of the energy in GeV. The resulting total number of CR, $N_{tot}$, is

$N_{tot}=N_{exi}+N_{imp}.$ (34)

Figure 12 reports the total number of CR as evaluated in Equation (34).

3.5. Variable Number of Injected Particles

We now analyse the case where the number of injected particles at $t=0$ varies as a function of the coefficient of diffusion, analysing the following 3D solution

$N\left(r\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}D\mathrm{,}\beta \mathrm{,}{D}_0\mathrm{,}{N}_0\right)=\frac{N_0{\left(\frac{D}{D_0}\right)}^{\beta }\sqrt{4}{\text{e}}^{-\frac{r^2}{4Dt}}}{16{\left(\pi Dt\right)}^{\frac{3}{2}}}\mathrm{,}$ (35)

Figure 9. The CR proton spectrum as measured by Voyager 1.

Figure 10. The CR proton spectrum as measured by Voyager 1 (empty stars) and polynomial fit of fourth degree.

where $N_0$ is the number of injected particles when $D=D_0$ and $\beta$ is an exponent which characterizes the transient diffusion. This equation, when $D=D_0$, becomes the well known 3D solution for the impulsive case, see Equation (16). The above number of diffusing particles has a maximum at

$D=-\frac{1}{2t\left(2\beta -3\right)},$ (36)

Figure 11. Number of particles in the case of 3D impulsive injection as function of energy at 1 au when $N=0.45\times {10}^{-11}$, $d=3$, $Z=1$, $B_{-6}=6000$ and $r=4.84\times {10}^{-6}\text{pc}$.

Figure 12. Number of existing CR (red full line) and the sum of existing CR and diffusing CR (blue dotted line).

which is $D=5.55\times {10}^{-4}$ with the parameters of Figure 13. According to Equation (27), the version in energy of solution (35) is

$N\left(r\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{E}_{\text{GeV}}\mathrm{,}\beta \mathrm{,}{E}_{\text{GeV},\text{0}}\mathrm{,}{N}_0\right)=\frac{8.67\times {10}^8{N}_0{\left(\frac{E_{\text{GeV}}}{E_{\text{GeV}\mathrm{,0}}}\right)}^{\beta }\sqrt{4}{\text{e}}^{-\frac{4.53\times {10}^6{r}^2{B}_{-6}}{E_{\text{GeV}}t}}}{{\left(\frac{E_{\text{GeV}}t}{B_{-6}}\right)}^{\frac{3}{2}}}\mathrm{,}$ (37)

Figure 13. The number of diffusing particles in the 3D impulsive case as function of a variable diffusion coefficient when $N_0=1$, $t=100$, $r=1$, $D_0=1$ and $\beta =-3$.

where $N_0$ is the number of injected particles when $E_{\text{GeV}}=E_{\text{GeV}\mathrm{,0}}$. The distribution in energy has a maximum at

$E=-\frac{9.068\times {10}^6{r}^2{B}_{-6}}{t\left(2\beta -3\right)}\text{GeV},$ (38)

and therefore equating the experimental energy where the energy has a maximum in the number of particles, $E_{\text{mode},\text{GeV}}$, and the theoretical one, allows introducing a relation between the involved parameters. As an example when $\beta =-3$, $Z=1$ and $r=4.848\times {10}^{-6}\text{pc}$ (1 au)

$t=\frac{2.3684\times {10}^{-5}{B}_{-6}}{E_{\exp }}\text{years}\mathrm{.}$ (39)

This relation allows eliminating the time from Equation (37), obtaining

$\begin{array}{l}N\left(r\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{E}_{\text{GeV}}\mathrm{,}\beta \mathrm{,}{E}_{\text{GeV}\mathrm{,0}}\mathrm{,}{E}_{\text{mode},\text{GeV}}\mathrm{,}{N}_0\right)\\ =\frac{8.67\times {10}^8{N}_0{\left(\frac{E_{\text{GeV}}}{E_{\text{GeV}\mathrm{,0}}}\right)}^{\beta }\sqrt{4}{\text{e}}^{\frac{0.5E_{\text{mode},\text{GeV}}\left(2\beta -3\right)}{E_{\text{GeV}}}}}{{\left(-\frac{9.068\times {10}^6{E}_{\text{GeV}}{r}^2}{E_{\text{mode},\text{GeV}}\left(2\beta -3\right)}\right)}^{\frac{3}{2}}}\mathrm{.}\end{array}$ (40)

The distribution in energy of this formula is reported in Figure 14 for the case of the Voyager 1 spacecraft and in Figure 15 for the whole spectrum of CR [54].

3.6. The Knee and the Second Knee for CR

The Kascade experiment [59] has measured the CR with energies in the range

Figure 14. The number of diffusing particles in the 3D impulsive case, see Equation (40) as function of the energy when $N_0=5\times {10}^{-12}$, $r=4.84\times {10}^{-6}\text{pc}$, $E_{\text{GeV}\mathrm{,0}}=2.6\times {10}^{-3}$, $E_{\text{mode},\text{GeV}}=6.2\times {10}^{-2}$, $\beta =-0.1$, (full red line) and Voyager 1 data, (empty green stars).

Figure 15. The number of diffusing particles, see Equation (40) as function of the energy when $N_0=1.0\times {10}^5$, $r=4.84\times {10}^{-6}\text{pc}$, $E_{\text{GeV}\mathrm{,0}}=0.21$, $E_{\text{mode},\text{GeV}}=0.21$, $\beta =-1$, (full red line) and CR distribution, (dotted green line).

[2.23 × 10⁶ GeV - 1.25 × 10⁹ GeV]. The data are available at Cosmic-Ray DataBase (CRDB), which has the address https://lpsc.in2p3.fr/crdb/; we selected the H-He-group, see Figure 16.

We now test the hypothesis that the CR in the range of energy of the Kascade

Figure 16. Flux of H-He versus energy per nucleus in Gev: experimental data (empty stars).

experiment are produced on the expanding surface of the local bubble [60] [61], here approximated by a sphere of radius $a=100\text{pc}$. Therefore the CR diffuses toward the inside of the sphere with a solution as reported in formula (23); we recall that the number of particles injected at a given diffusion coefficient is assumed to be constant over time. Conversely the number of particles, the boundary concentration at $r=a$, as a function of the energy, is assumed to scale as

$f=f_a{\left(\frac{E_{\min }}{E}\right)}^{\alpha },$ (41)

where $f_a$ is the concentration at $\left(\mathrm{0,}a\right)$ and $E_{\min }$ the minimum energy considered. We now assume that for a preexisting spectrum of energy, $N_{pre}\left(E\right)$, two events of diffusion for different radii of the ball, $N_{e\mathrm{,1}}\left(E\right)$ and $N_{e\mathrm{,2}}\left(E\right)$, are summed in order to obtain the total number of CR, $N_{tot}$,

$N_{tot}\left(E\right)=N_{pre}\left(E\right)+N_{e,1}\left(E\right)+N_{e,2}\left(E\right).$ (42)

The preexisting number of CR is assumed to be $N_{pre}=745203.5\times E^{-2.97}$ and the parameters of the two events of diffusion are reported in Table 1. Figure 17 reports $N_{tot}$ as well the data of the Kascade experiment.

3.7. The ankle for CR

The ankle in the distribution of high-energy CR is at ≈10⁷TeV and characterizes the transition from galactic to extra-galactic CR [62]; Figure 18 reports the experimental energy distribution according to Figure 1 in [63].

We now test solution (40), which represents an impulsive 3D solution in the presence of a variable number of injected particles situated at $r=20\text{pc}$, which

Table 1. Two events of diffusion inside the local bubble when $u_0=f/1000$, $r=20\text{pc}$, $B_{-6}=1.8$, $Z=1$, $d=3$ and the number of iterations of the series is 20.

Figure 17. $N_{tot}\left(E\right)$ as function of energy (full line) with data as in Table 1 and data of the Kascade experiment (empty stars).

Figure 18. CR energy in TeV: experimental data (empty stars) according to [63].

is the distance of the expanding surface of the local bubble. Figure 19 reports the existing number of CR, $N_{exi}$, evaluated as in Equation (40). We now add to an existing power law distribution an impulsive phenomena; Figure 20 reports the number of CR, $N_{imp}$, in a burst at $r=20\text{pc}$.

The total number of CR as evaluated with formula (34) is reported in Figure 21, which thus makes visible a theoretical explanation for the ankle.

Figure 19. Flux of H versus energy per nucleus in Gev: experimental data (empty green stars) according to [63] and theoretical existing distribution (red full line) when $N_0=2\times {10}^{-7}$, $d=3$, $Z=1$, $E_{\text{mode},\text{GeV}}=10888.59$, $\beta =-1.5$, and $r=20\text{pc}$.

Figure 20. Number of particles in a burst for a 3D impulsive injection as function of the energy when $N=3.0\times {10}^{-16}$, $d=3$, $Z=1$, $B_{-6}=6000$, $E_{\text{mode},\text{GeV}}={10}^{10}$, $\beta =-1.5$, and $r=20\text{pc}$.

Figure 21. Number of existing CR (green points) and the sum of existing CR and diffusing CR (red line).

4. Conclusions

PDE & Boundary Conditions

Two new solutions for the diffusion equation have been derived: the first one was derived for the case of a 1D diffusion and an exponential profile for the number of particles, see Equation (19); the second one gives the 2D diffusion for the case of a power law profile for the number of particles, see Equation (22).

CR and Diffusion

In order to obtain a distribution in energy for cosmic rays (CR) in 3D which can be compared to the observed one, the number of injected particles should be a function of the diffusion coefficient, see Equation (35), or its energy equivalent, see Equation (37). The theoretical distribution has a maximum at the value of the diffusion coefficient given by Equation (36) or at the value of energy represented by Equation (38).

Knee and Ankle

In the framework of summing an existing distribution in energy for the CR originating from the local bubble with that originating from other sources, see Equation (42), it is possible to simulate the knee and the second knee, see Figure 17. The ankle is simulated in Figure 19 through the superposition of two events of solar origin assuming a time independent solution for the energy as given by Equation (40).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Hess, V. (2018) On the Observations of the Penetrating Radiation during Seven Balloon Flights.
[2]	Andrews, M.R. (1925) High Frequency Rays of Cosmic Origin. Science, 62, 445-448. https://doi.org/10.1126/science.62.1612.445
[3]	Wright, C.S. (1926) Cosmic Rays. Nature, 117, 54-56. https://doi.org/10.1038/117054a0
[4]	Millikan, R.A. (1926) High Frequency Rays of Cosmic Origin. Popular Astronomy, 34, 232-338. https://doi.org/10.1073/pnas.12.1.48
[5]	Millikan, R.A. (1926) Cosmic Rays. Scientific American, 134, 149. https://doi.org/10.1038/scientificamerican0326-149
[6]	Parker, E.N. (1958) Cosmic-Ray Modulation by Solar Wind. Physical Review, 110, 1445-1449. https://doi.org/10.1103/PhysRev.110.1445
[7]	Ehmert, A. (1960) On the Modulation of Primary Cosmic Ray Spectrum by Solar Activity. International Cosmic Ray Conference, Vol. 4, 142.
[8]	Dorman, L.I. (1960) On the Theory of the Modulation of Cosmic Rays by the Solar Wind. International Cosmic Ray Conference, Vol. 4, 320.
[9]	Potgieter, M.S. (2013) Solar Modulation of Cosmic Rays. Living Reviews in Solar Physics, 10, Article No. 3. https://doi.org/10.12942/lrsp-2013-3
[10]	Miyake, S., Kataoka, R. and Sato, T. (2017) Cosmic Ray Modulation and Radiation Dose of Aircrews during the Solar Cycle 24/25. Space Weather, 15, 589-605. https://doi.org/10.1002/2016SW001588
[11]	Neher, H.V. (1961) Cosmic-Ray Knee in 1958. Journal of Geophysical Research, 66, 4007. https://doi.org/10.1029/JZ066i012p04007
[12]	Hillas, A.M. (1979) The Knee of the Cosmic-Ray Spectrum: Not a Magnetic Trapping Effect? International Cosmic Ray Conference, Vol. 8, 7.
[13]	Forman, M.A. (1981) Could the Cosmic-Ray Knee at 1015 eV Be Due to the Finite Time for Supernova Shock Acceleration? Bulletin of the American Astronomical Society, 13, 796.
[14]	Clay, R.W. (1984) Cosmic Ray Anisotropy and the Knee of the Energy Spectrum. Australian Journal of Physics, 37, 97-103. https://doi.org/10.1071/PH840097
[15]	Clay, R.W. (1988) The Knee of the Cosmic Ray Energy Spectrum. Australian Journal of Physics, 41, 729-733. https://doi.org/10.1071/PH880729
[16]	Kang, H., Rachen, J.P. and Biermann, P.L. (1996) Contributions to the Cosmic Ray Flux above the Ankle: Clusters of Galaxies. Journal of Korean Astronomical Society Supplement, 29, S271-S272.
[17]	Kang, H., Rachen, J.P. and Biermann, P.L. (1997) Contributions to the Cosmic Ray Flux above the Ankle: Clusters of Galaxies. MNRAS, 286, 257-267. https://doi.org/10.1093/mnras/286.2.257
[18]	Gelmini, G. and Varieschi, G. (2002) Cosmic Rays above the Ankle from Z-Bursts.
[19]	Compton, A.H. and Getting, I.A. (1935) An Apparent Effect of Galactic Rotation on the Intensity of Cosmic Rays. Physical Review, 47, 817-821. https://doi.org/10.1103/PhysRev.47.817
[20]	Vallarta, M.S. and Feynman, R.P. (1939) The Scattering of Cosmic Rays by the Stars of a Galaxy. Physical Review, 55, 506-507. https://doi.org/10.1103/PhysRev.55.506.2
[21]	Kiepenheuer, K.O. (1950) Cosmic Rays as the Source of General Galactic Radio Emission. Physical Review, 79, 738-739. https://doi.org/10.1103/PhysRev.79.738
[22]	Unsöld, A. (1951) Cosmic Radiation and Cosmic Magnetic Fields. I. Origin and Propagation of Cosmic Rays in Our Galaxy. Physical Review, 82, 857-863. https://doi.org/10.1103/PhysRev.82.857
[23]	Burbidge, G.R. (1957) Acceleration of Cosmic-Ray Particles among Extragalactic Nebulae. Physical Review, 101, 269-271. https://doi.org/10.1103/PhysRev.107.269
[24]	Fujimoto, Y., Hasegawa, H. and Taketani, M. (1964) Part I. Cosmic Rays in Galactic and Extragalactic Space. Progress of Theoretical Physics Supplement, 30, 32-85. https://doi.org/10.1143/PTPS.30.32
[25]	Laster, H. (1964) Galactic and Extragalactic Propagation of Cosmic Rays. Physical Review, 135, 1274-1279. https://doi.org/10.1103/PhysRev.135.B1274
[26]	Fermi, E. (1949) On the Origin of the Cosmic Radiation. Physical Review, 75, 1169-1174. https://doi.org/10.1103/PhysRev.75.1169
[27]	Bohm, D. and Gross, E.P. (1948) Plasma Oscillations as a Cause of Acceleration of Cosmic-Ray Particles. Physical Review, 74, 624. https://doi.org/10.1103/PhysRev.74.624
[28]	Parker, E.N. (1955) Hydromagnetic Waves and the Acceleration of Cosmic Rays. Physical Review, 99, 241-263. https://doi.org/10.1103/PhysRev.99.241
[29]	Thompson, W.B. (1955) On the Acceleration of Cosmic-Ray Particles by Magneto-Hydrodynamic Waves. Proceedings of the Royal Society of London Series A, 233, 402-406. https://doi.org/10.1098/rspa.1955.0275
[30]	Laster, H. (1957) Factors Affecting the Galactic Diffusion of Cosmic Rays. Physical Review, 107, 1112. https://doi.org/10.1103/PhysRev.107.1112
[31]	Jokipii, J.R. (1967) Cosmic-Ray Propagation. II. Diffusion in the Interplanetary Magnetic Field. ApJ, 149, 405. https://doi.org/10.1086/149265
[32]	Burlaga, L.F. (1967) Anisotropic Diffusion of Solar Cosmic Rays. Journal of Geophysical Research, 72, 4449-4466. https://doi.org/10.1029/JZ072i017p04449
[33]	Jokipii, J.R. (1968) Backscatter and Diffusion of Solar Cosmic Rays. The Astronomical Journal Supplement, 73, 66. https://doi.org/10.1086/149612
[34]	Miroshnichenko, L.I. (1969) Solar Cosmic-Ray Diffusion in the Presence of a Magnetic Barrier in Interplanetary Space. Geomagnetism and Aeronomy, 9, 271.
[35]	Bel’Skii, S.A. (1970) The Magnetic Fields of Meteor Streams as a Possible Mechanism for the Diffusion of Cosmic Rays in Circumsolar Space. Astronomicheskii Zhurnal, 47, 201-205.
[36]	Jokipii, J.R. and Parker, E.N. (1970) On the Convection, Diffusion, and Adiabatic Deceleration of Cosmic Rays in the Solar Wind. ApJ, 160, 735. https://doi.org/10.1086/150465
[37]	Lingenfelter, R.E., Ramaty, R. and Fisk, L.A. (1971) Compound Diffusion of Cosmic Rays. Astrophysical Letters, 8, 93-97.
[38]	Jones, F.C., Birmingham, T.J. and Kaiser, T.B. (1973) Investigation of Resonance Integrals Occurring in Cosmic-Ray Diffusion Theory. ApJ Letters, 180, L139. https://doi.org/10.1086/181170
[39]	Bagdasarian, M.B. and Dorman, L.I. (1980) Cosmic Ray Modulation by a Shock Wave Taking into Account Adiabatic Cooling and Transversal Diffusion (Automodel Solution). Geomagnetism and Aeronomy, 20, 990.
[40]	Bazilevskaya, G.A. and Golynskaya, R.M. (1990) Features of the Focused Diffusion of Solar Cosmic Rays. Geomagnetism and Aeronomy, 30, 725-727.
[41]	Candia, J., Roulet, E. and Epele, L.N. (2002) Turbulent Diffusion and Drift in Galactic Magnetic Fields and the Explanation of the Knee in the Cosmic Ray Spectrum. Journal of High Energy Physics, 2002, 33. https://doi.org/10.1088/1126-6708/2002/12/033
[42]	Shalchi, A., Büsching, I., Lazarian, A. and Schlickeiser, R. (2010) Perpendicular Diffusion of Cosmic Rays for a Goldreich—Sridhar Spectrum. ApJ, 725, 2117-2127. https://doi.org/10.1088/0004-637X/725/2/2117
[43]	Medvedev, M.V. and Medvedev, V.V. (2015) Asymmetric Diffusion of Cosmic Rays. Physics of Plasmas, 22, Article ID: 091504. https://doi.org/10.1063/1.4928942
[44]	Wiengarten, T., Oughton, S., Engelbrecht, N.E., Fichtner, H., Kleimann, J. and Scherer, K. (2016) A Generalized Two-Component Model of Solar Wind Turbulence and ab Initio Diffusion Mean-Free Paths and Drift Lengthscales of Cosmic Rays. ApJ, 833, 17. https://doi.org/10.3847/0004-637X/833/1/17
[45]	Buonocore, S. and Sen, M. (2021) Anomalous Diffusion of Cosmic Rays: A Geometric Approach. AIP Advances, 11, Article ID: 055221. https://doi.org/10.1063/5.0049401
[46]	Gould, H. and Tobochnik, J. (1988) An Introduction to Computer Simulation Methods. Addison-Wesley, Reading.
[47]	Olver, F.W.J., Lozier, D.W., Boisvert, R.F. and Clark, C.W. (2010) NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge.
[48]	Lang, K.R. (1999) Astrophysical Formulae. Third Edition, Springer, New York. https://doi.org/10.1007/978-3-662-21639-2
[49]	Longair, M.S. (2011) High Energy Astrophysics III ed. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511778346
[50]	Hillas, A.M. (1984) The Origin of Ultra-High-Energy Cosmic Rays. ARA&A, 22, 425. https://doi.org/10.1146/annurev.aa.22.090184.002233
[51]	Adriani, O., Barbarino, G.C., Bazilevskaya, G.A., et al. (2011) PAMELA Measurements of Cosmic-Ray Proton and Helium Spectra. Science, 332, 69-72. https://doi.org/10.1126/science.1199172
[52]	Aguilar, M., Aisa, D., Alpat, B., et al. (2015) Precision Measurement of the Proton Flux in Primary Cosmic Rays from Rigidity 1 GV to 1.8 TV with the Alpha Magnetic Spectrometer on the International Space Station. Physical Review Letters, 114, Article ID: 171103.
[53]	Aguilar, M., Aisa, D., Alpat, B., et al. (2015) Precision Measurement of the Helium Flux in Primary Cosmic Rays of Rigidities 1.9 GV to 3 TV with the Alpha Magnetic Spectrometer on the International Space Station. Physical Review Letters, 115, Article ID: 211101.
[54]	Zyla, P., et al. (Particle Data Group) (2020) Review of Particle Physics. PTEP, 2020, 083C01. https://doi.org/10.1093/ptep/ptaa104
[55]	Stone, E.C., Cummings, A.C., McDonald, F.B., Heikkila, B.C., Lal, N. and Webber, W.R. (2013) Voyager 1 Observes Low-Energy Galactic Cosmic Rays in a Region Depleted of Heliospheric Ions. Science, 341, 150-153. https://doi.org/10.1126/science.1236408
[56]	Potgieter, M. (2014) Very Local Interstellar Spectra for Galactic Electrons, Protons and Helium. Brazilian Journal of Physics, 44, 581-588. https://doi.org/10.1007/s13538-014-0238-2
[57]	Cholis, I., Hooper, D. and Linden, T. (2016) A Predictive Analytic Model for the Solar Modulation of Cosmic Rays. Physical Review D, 93, Article ID: 043016. https://doi.org/10.1103/PhysRevD.93.043016
[58]	Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (1992) Numerical Recipes in Fortran. The Art of Scientific Computing. Cambridge University Press, Cambridge.
[59]	Schoo, S., Apel, W.D., Arteaga-Velázquez, J.C., et al. (2015) The Energy Spectrum of Cosmic Rays in the Range from 10{14} to 10{18} eV. 34th International Cosmic Ray Conference, Vol. 34, 263.
[60]	Zaninetti, L. (2020) On the Shape of the Local Bubble. International Journal of Astronomy and Astrophysics, 10, 11-27. https://doi.org/10.4236/ijaa.2020.101002
[61]	Zaninetti, L. (2021) Energy Conservation in the Thin Layer Approximation: VI. Bubbles and Super-Bubbles. International Journal of Astronomy and Astrophysics, 11, 370-391. https://doi.org/10.4236/ijaa.2021.113017
[62]	Apel, W.D., Arteaga-Velàzquez, J.C., Bekk, K., et al. (2013) Ankle-Like Feature in the Energy Spectrum of Light Elements of Cosmic Rays Observed with KASCADE-Grande. Physical Review D, 87, Article ID: 081101. https://doi.org/10.1103/PhysRevD.87.081101
[63]	Schroeder, F., AbuZayyad, T., Anchordoqui, L., et al. (2019) High-Energy Galactic Cosmic Rays. Bulletin of the American Astronomical Society, 51, 131.