Transport in Astrophysics: I. Diffusion of Solar and Galactic Cosmic Rays ()
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^{7} 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^{18} eV and 10^{20} 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,
$\stackrel{\xaf}{{R}^{2}\left(t\right)}$, according to Equation (8.38) in [46] is
$\stackrel{\xaf}{{R}^{2}\left(t\right)}=2dDt\text{\hspace{1em}}\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{{\stackrel{\xaf}{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{\text{\hspace{0.17em}}}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{\hspace{0.05em}}{\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}_{\mathrm{max}}$ where
${t}_{\mathrm{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\text{\hspace{0.17em}}}\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\text{\hspace{0.17em}}}\frac{{\text{pc}}^{\text{2}}}{\text{yr}}$ (green line), and for the error solution solution when
${N}_{0}=1$,
$D=\mathrm{1\text{\hspace{0.17em}}}\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}_{\mathrm{max}}$ where
${t}_{\mathrm{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{;}D\right)=\frac{{N}_{0}\sqrt{4}\text{\hspace{0.05em}}{\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}_{\mathrm{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\mathrm{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\text{\hspace{0.17em}}}\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\mathrm{sin}\left(\frac{n\pi x}{L}\right)\mathrm{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)}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}.\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{\hspace{0.05em}}\text{\hspace{0.05em}}\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(2{N}_{0}-2{N}_{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}}}\mathrm{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{\hspace{0.17em}}\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.081{E}_{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^{15} eV, and
${B}_{-6}$ is the magnetic field expressed in 10^{−6} 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.055134{E}_{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^{5} GeV]:
${I}_{H}\left(E\right)\approx 1.24\mp {10}^{4}{\left(\frac{E}{1\text{\hspace{0.17em}}\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}_{\mathrm{min}}$ and
${D}_{\mathrm{max}}$ is
$\frac{N\left(r,t;{D}_{\mathrm{min}}\right)}{N\left(r,t;{D}_{\mathrm{max}}\right)}=\frac{{\text{e}}^{-\frac{{r}^{2}\left({D}_{\mathrm{max}}-{D}_{\mathrm{min}}\right)}{4{D}_{\mathrm{min}}t{D}_{\mathrm{max}}}}{D}_{\mathrm{max}}^{\frac{3}{2}}}{{D}_{\mathrm{min}}^{\frac{3}{2}}}.$ (29)
This ratio can be parametrized as
$\frac{N\left(r,t;{D}_{\mathrm{min}}\right)}{N\left(r,t;{D}_{\mathrm{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}_{\mathrm{min}}=0.02$,
${D}_{\mathrm{max}}=0.2$ and
${N}_{0}=1$ is
$t=1.0857{r}^{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 ={\mathrm{log}}_{10}\left(\frac{N\left(r,t;{D}_{\mathrm{min}}\right)}{N\left(r,t;{D}_{\mathrm{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{\hspace{0.17em}}\text{GeV}}\right)}^{-3}$ when
${D}_{\mathrm{min}}=0.02$,
${D}_{\mathrm{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{\hspace{0.17em}}\text{GeV}}\right)}^{\alpha}$, when
$r=10\text{\hspace{0.17em}}\text{pc}$,
${D}_{\mathrm{min}}=0.02$ and
${D}_{\mathrm{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^{−3} 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{;}D\mathrm{,}\beta \mathrm{,}{D}_{0}\mathrm{,}{N}_{0}\right)=\frac{{N}_{0}{\left(\frac{D}{{D}_{0}}\right)}^{\beta}\sqrt{4}\text{\hspace{0.05em}}{\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{\hspace{0.05em}}\text{\hspace{0.05em}}\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{;}{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{\hspace{0.05em}}{\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{\hspace{0.05em}}\text{pc}$ (1 au)
$t=\frac{2.3684\times {10}^{-5}{B}_{-6}}{{E}_{\mathrm{exp}}}\text{years}\mathrm{.}$ (39)
This relation allows eliminating the time from Equation (37), obtaining
$\begin{array}{l}N\left(r\mathrm{,}t\mathrm{;}{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{\hspace{0.05em}}{\text{e}}^{\frac{0.5{E}_{\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{\hspace{0.05em}}\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{\hspace{0.05em}}\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^{6} GeV - 1.25 × 10^{9} 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{\hspace{0.17em}}\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}_{\mathrm{min}}}{E}\right)}^{\alpha},$ (41)
where
${f}_{a}$ is the concentration at
$\left(\mathrm{0,}a\right)$ and
${E}_{\mathrm{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^{7}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{\hspace{0.17em}}\text{pc}$, which
Table 1. Two events of diffusion inside the local bubble when
${u}_{0}=f/1000$,
$r=20\text{\hspace{0.17em}}\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{\hspace{0.17em}}\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{\hspace{0.17em}}\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{\hspace{0.17em}}\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).