Particle Spectra and Mass Composition in the Ultra-High Energy Region in the Framework of the Galactic Origin of Cosmic Rays

TL;DR

This paper explores a Galactic origin model for ultra-high energy cosmic rays, accounting for source acceleration limits, inhomogeneous galactic matter, and magnetic fields, to explain observed spectra and composition variations.

Contribution

It presents a self-consistent Galactic origin scenario incorporating fractal-like matter distribution and Levy flights, aligning with recent experimental data on cosmic ray composition.

Findings

Model explains cosmic ray spectra up to 10^{20} eV.

Predictions align with recent composition measurements.

Scenario suggests Galactic sources can account for ultra-high energies.

Abstract

The possibility for a self-consistent description of all the basic features of the observed cosmic ray spectra and primary composition variations in the energy range of $10^{15}\div 10^{20}$ eV within the Galactic origin scenario is examined. We assume the existence of Galactic sources that accelerate particles up to $\sim 3\cdot 10^{18}Z$ eV and take into account a highly inhomogeneous (fractal-like) distribution of matter and magnetic fields in the Galaxy that leads to extremely large free paths of particles ("L\'{e}vy flights"), along with an overwhelming contribution to the cosmic ray fluxes observed above $\sim 10^{18}$ eV from particles reaching the Solar System without scattering. Our scenario was refined on the basis of recent experimental results on primary mass composition. Model predictions, which could be verified with the improved high-precision measurements in the nearest future are discussed.

Tables

Table 1. Table 1: Anomalous diffusion model parameters

Parameter	Value
$\alpha$	$1.7$
$\beta$	$0.8$
$D_0(\alpha ,\beta )$	$1.5\cdot {10}^{-3}$ pc1.7y-0.8
$\delta$	$0.27$
${\delta }_L$	$\delta /2$
$\gamma$	$2.85$
$E_0$	$4\cdot {10}^{18}Z$ eV
$T$	${10}^4$ y

Table 2. Table 2: The knee energies for different elemental groups of nuclei

Nuclei	Knee energy, eV
H	$\left(6.1÷6.5\right)\cdot {10}^{14}$
He	$\left(1.8÷2.2\right)\cdot {10}^{15}$
CNO	$\left(5.8÷6.3\right)\cdot {10}^{15}$
NeMgSi	$\left(0.8÷1.2\right)\cdot {10}^{16}$
Fe	$\left(2.6÷3.0\right)\cdot {10}^{16}$

Equations

$$J(\mathbf{r},t,E)=J_G(\mathbf{r},E)+J_L(\mathbf{r},t,E)+J_{NS}(\mathbf{r},E).$$

$$\begin{gathered} \text{leavevmode}\frac{\partial N(\mathbf{r},t,E)}{\partial t}=-D(E,\alpha ,\beta ){\mathrm{D}}_{0+}^{1-\beta }(-\Delta {)}^{\alpha /2}N(\mathbf{r},t,E)+ \\ +S(\mathbf{r},t,E). \end{gathered}$$

$$\begin{gathered} N(\mathbf{r},t,E)=\frac{S_0{E}^{-\gamma }}{D(E,\alpha ,\beta {)}^{3/\alpha }}\times \\ \underset{\max [0,t-T]}{\overset{t}{\int }}d\tau {\tau }^{-3\beta /\alpha }{\Psi }_3^{(\alpha ,\beta )}\left(|\mathbf{r}|(D(E,\alpha ,\beta ){\tau }^{\beta }{)}^{-1/\alpha }\right), \end{gathered}$$

$${\Psi }_3^{(\alpha ,\beta )}(\rho )=\underset{0}{\overset{\infty }{\int }}{g}_3^{(\alpha )}(r{\tau }^{\beta })q_1^{(\beta ,1)}(\tau ){\tau }^{3\beta /\alpha }d\tau .\text{vspace}\ast -3mm$$

$$\eta =-\frac{d\log N}{d\log E}=\gamma +\frac{\delta }{\alpha }\Xi ,$$

$$\begin{gathered} \Xi =3-\frac{2\pi r^2}{D(E,\alpha ,\beta {)}^{2/\alpha }}\times \\ \times \frac{\underset{\max [0,t-T]}{\overset{t}{\int }}d\tau {\tau }^{-5\beta /\alpha }{\Psi }_5^{(\alpha ,\beta )}\left(|\mathbf{r}|(D(E,\alpha ,\beta ){\tau }^{\beta }{)}^{-1/\alpha }\right)}{\underset{\max [0,t-T]}{\overset{t}{\int }}d\tau {\tau }^{-3\beta /\alpha }{\Psi }_3^{(\alpha ,\beta )}\left(|\mathbf{r}|(D(E,\alpha ,\beta ){\tau }^{\beta }{)}^{-1/\alpha }\right)}. \end{gathered}$$

$$\delta =\left(\eta {|}_{E\gg E_k}-\eta {|}_{E\ll E_k}\right)\frac{\beta }{1+\beta },\gamma =\eta {|}_{E\ll E_k}+\delta .$$

$$\gamma \sim 2.85÷2.95,\beta \sim \mathrm{0.8.}$$

$$J(\mathbf{r},t,E)=\frac{v}{4\pi}\Biggl{[}S_{G}E^{-\gamma-\delta/\beta}+\frac{S_{0}E^{-\gamma}}{D(E,\alpha,\beta)^{3/\alpha}}\times\Biggr{.}\\ \times\sum\limits_{\begin{subarray}{c}r_{j}<1\;\text{kpc}\\ t_{j}<10^{6}\;\text{yr}\end{subarray}}\int\limits_{\max[0,t_{j}-T]}^{t_{j}}d\tau\tau^{-3\beta/\alpha}\Psi_{3}^{(\alpha,\beta)}\left(|\mathbf{r}_{j}|(D(E,\alpha,\beta)\tau^{\beta})^{-1/\alpha}\right)+\\[-14.22636pt] +\Biggl{.}S_{NS}\sum\limits_{r_{j}<1\;\text{kpc}}E^{-\gamma+\delta_{L}}|\mathbf{r}_{j}|^{-\alpha}\Biggr{]}\exp\left(-\frac{E}{E_{0}}\right).$$

Full text

We assume the existence of Galactic sources that accelerate particles up to $\sim 3\cdot 10^{18}Z$ eV and take into account a highly inhomogeneous (fractal-like) distribution of matter and magnetic fields in the Galaxy that leads to extremely large free paths of particles (“Lévy flights”), along with an overwhelming contribution to the cosmic ray fluxes observed above $\sim 10^{18}$ eV from particles reaching the Solar System without scattering. Our scenario was refined on the basis of recent experimental results on primary mass composition. Model predictions, which could be verified with the improved high-precision measurements in the nearest future are discussed.

1 Introduction

Physical mechanisms governing the features of energy spectra of cosmic rays (CRs) and their composition variations in the energy range that spans from PeV to tens of EeV are still under debate. The so-called standard scenario assumes that the Galactic component of CRs accelerated by supernova remnants extends up to energies of a few $10^{17}$ eV with gradually heavier composition due to the rigidity-dependent maximum acceleration energy of nuclei. Then Galactic CRs merge into a light and flatter extragalactic component, which terminates at $\sim 5\cdot 10^{19}$ eV due to the GZK effect.

Recent measurements by several new generation experiments and also the results of some new methods of experimental data analysis Apel-a:2013 ; Buitink:2016 ; Settimo:2016 ; Aab:2016 ; Abbasi:2015 have severely challenged the standard scenario. In particular, an ankle-like feature in the light CR component measured by KASCADE-Grande at $E>10^{17}$ eV Apel-a:2013 together with evidence for a light-mass fraction (protons and helium) of about 80 percent at $10^{17}\div 10^{17.5}$ eV Buitink:2016 cause difficulties for an explanation in the framework of the standard model without strong additional assumptions. At ultra-high energies the chemical composition is still ambiguous. Mixed composition including nuclei heavier than helium is reported by the Pierre Auger Observatory (PAO) on the basis of various improved methods and techniques Settimo:2016 ; Aab:2016 , while the Telescope Array (TA) interpretes their data as being compatible with a largely protonic composition Abbasi:2015 . Thus, the non-GZK origin of the energy spectrum suppression at the highest energies (e.g. the acceleration limits in Galactic cosmic ray sources above $10^{18}$ eV) can not be excluded.

In our papers Lagutin:2005ijmp ; Lagutin:2015 ; Volkov:2015 it was shown that basic structures in the all particle CR spectrum (knee, ankle and cut-off) can be reproduced under the assumption that essentially all CRs, including those at the highest energies, originate from Galactic sources, which can accelerate particles up to the cut-off energy. The highly inhomogeneous (fractal-like) distribution of matter and magnetic fields in the Galaxy that leads to extremely large free paths of particles (so called “Lévy flights”), along with large contribution of non-scattered particles to cosmic ray fluxes above $\sim 10^{18}$ eV is the main element of our model – the anomalous diffusion model. At the present stage the new more precise data, mainly the spectrum features measured separately for light and heavy components along with the improved composition estimations, provide an opportunity for refinements and updates of the model.

In the present paper we revise the possibility of the self-consistent description of the CRs energy spectra and composition results within the Galactic cosmic ray origin scenario based on the anomalous diffusion model and discuss the crucial model predictions, which could be verified with improved high-precision measurements in the near future.

2 Key elements of the proposed scenario

The following basic principles and assumptions were initially included in the proposed scenario.

The existence of Galactic sources that accelerate particles with atomic charge $Z$ up to $E_{0}\sim 4\cdot 10^{18}Z$ eV. Consequently, the observed suppression of the primary CR spectrum at $\sim 5\cdot 10^{19}$ eV is not an extragalactic GZK feature; it reflects the acceleration limit of Galactic CR sources. 2. 2.

The highly inhomogeneous (fractal-like) distribution of matter and magnetic fields in the Galaxy that leads to the anomalous diffusion of CRs manifested, in particular, by abnormally large free paths of particles (so-called “Lévy flights”) with a power-law distribution. 3. 3.

An overwhelming contribution to cosmic ray fluxes observed above $\sim 10^{18}$ eV from particles reaching the Solar system without scattering.

Note, that the main argument contrary to any model assuming the Galactic origin of ultra-high energy cosmic rays is that particles originating from Galactic sources could hardly be as isotropic as presently observed. However, in recent works Kumar:2014 ; Eichler:2016 it was shown that the isotropy can be understood in the context of Galactic production if a sufficiently careful treatment of CR propagation in the Galaxy is undertaken.

The spatial distribution of sources also suggests the separation of the observed CR flux into three components as follows:

[TABLE]

Here

• •

$J_{G}$ is the global spectrum component determined by the multiple old ($t\geq 10^{6}$ yr) distant ($r\geq 1$ kpc) sources.

• •

$J_{L}$ is the local component, i.e. the contribution of nearby ($r<1$ kpc) young ($t<10^{6}$ yr) sources. The spatial and temporal coordinates of the local sources are presented in LagutinTymentcev:2004 .

• •

$J_{NS}$ is the flux of non-scattered particles.

We assume that the non-scattered component is also formed by nearby ($r<1$ kpc) sources, defining the spectrum in the ultra-high energy region, and provides the observed flattening of the spectrum at $E\geq 10^{18}$ eV.

The flux of non-scattered particles, $J_{NS}$, is determined by the injected flux $(S_{0}E^{-\gamma}\exp(-E/E_{0})$, where $E_{0}=4\cdot 10^{18}Z$ eV), multiplied by a factor describing the probability of reaching the observer at a distance, $r$, from the source for a particle with a certain energy and atomic number without scattering. This factor has a power-law asymptotic with respect to $r$ (the Lévy flight probability; see the next section for the explicit form) and also suppressed at energies $E<3\cdot 10^{17}Z$ eV reflecting the fact that for the nucleus with gyroradius less than the typical Galactic inhomogeneity size this inhomogeneity would be opaque.

3 Parameters of the cosmic ray anomalous diffusion model

The highly inhomogeneous character of matter distribution and associated magnetic fields in the Galaxy should be adequately incorporated into the cosmic ray diffusion model. A physically reasonable way for the generalization of the normal diffusion model is to abandon the assumption about statistical homogeneity of the distribution of matter in favour of its fractal distribution. A principle consequence of this generalization is the power-law distribution of free paths $r$ in such a medium $p(\mathbf{r},E)\propto A(E,\alpha)r^{-\alpha-1},r\rightarrow\infty,0<\alpha<2$ — so-called Lévy flights. Besides, the intermittent magnetic field of the fractal-like interstellar medium (ISM) leads to a higher probability of a long stay of particles in inhomogeneities, leading to a presence of the so-called Lévy traps. In the general case, the probability density function $q(t,E)$ of time $t$, during which a particle is trapped in the inhomogeneity (Lévy trap), also has a power-law behaviour: $q(t,E)\propto B(E,\beta)t^{-\beta-1},t\rightarrow\infty,\beta<1$.

Generalization of the homogeneous normal diffusion model to the case of inhomogeneous (fractal-like) ISM, has been made for the first time in our papers Lagutin:2000 ; Lagutin:2001 . Later, it was shown Lagutin:2003 ; Lagutin:2003ew ; LagutinTymentcev:2004 ; Lagutin:2009 that an anomalous cosmic ray diffusion model allows to describe the main features of nuclei, electron and positron spectra observed in the Solar system. Particularly, in the anomalous diffusion model the key feature of the all particle energy spectrum — the knee at $3\cdot 10^{15}$ eV — appears naturally without additional assumptions.

The equation for the density of particles with energy $E$ at the location $\mathbf{r}$ and time $t$, generated in a fractal-like medium by Galactic sources with a distribution density $S(\mathbf{r},t,E)$ can be written as Lagutin:2003 ; LagutinTymentcev:2004

[TABLE]

Here $\mathrm{D}_{0+}^{1-\beta}$ denotes the Riemann-Liouville fractional derivative Samko:1993 and $(-\Delta)^{\alpha/2}$ is the fractional Laplacian (“Riesz operator”) Samko:1993 . The anomalous diffusion coefficient $D(E,\alpha,\beta)\sim A(E,\alpha)/B(E,\beta)=D_{0}(\alpha,\beta)E^{\delta}$.

The solution of Eq. (1) for a point impulse source with a power-law injection spectrum and emission time $T$ $S(\mathbf{r},t,E)=S_{0}E^{-\gamma}\delta(\mathbf{r})\Theta(T-t)\Theta(t)$ ($\Theta(\tau)$ is the step function) has the form Lagutin:2003 ; LagutinTymentcev:2004

[TABLE]

where $\Psi_{3}^{(\alpha,\,\beta)}(\rho)$ is the density of the fractional stable distribution Uchaikin:1999a ; Zolotarev:1999

[TABLE]

Using the representation $N=N_{0}E^{-\eta}$ and the property $d\Psi_{m}^{(\alpha,\,\beta)}(\rho)/d\rho=-2\pi\rho\Psi_{m+2}^{(\alpha,\,\beta)}(\rho)$ of the scaling function Uchaikin:1999a , one can easy find the spectral exponent $\eta$ for observed particles:

[TABLE]

where

[TABLE]

Let $E_{k}$ be a solution of the equation $\Xi(E)=0$. One can see from (3) and (4) that at $E=E_{k}$ the spectral exponent for observed particles $\eta$ is equal to the spectral exponent for particles generated by the source: $\eta(E_{k})=\gamma$. Since the exponent $\eta|_{E\ll E_{k}}=\gamma-\delta$ is less than $\gamma$ at $E\ll E_{k}$, but the exponent $\eta|_{E\gg E_{k}}=\gamma+\delta/\beta$ is greater than $\gamma$ at $E\gg E_{k}$, $E_{k}$ can be called the “knee” energy.

From experimental values of $\eta|_{E\ll E_{k}}$ and $\eta|_{E\gg E_{k}}$ one can derive the main parameters of the model $(\gamma,\delta)$ versus the spectral exponent $\beta$ of the “Lévy waiting time”:

[TABLE]

Since $\eta|_{E\gg E_{k}}-\eta|_{E\ll E_{k}}\sim 0.6$, $\eta|_{E\ll E_{k}}\sim 2.55\div 2.65$ Bartoli:2015 , $\delta\sim 0.27$, the last equations permit to retrieve self-consistently both spectral exponents $\gamma$ and $\beta$:

[TABLE]

We note that a similar steep spectrum of accelerated particles has been observed from supernova remnants W44 Abdo:2010 and IC 443 Abdo:2010_1 with the Fermi Large Area Telescope, and W49B with H.E.S.S. and Fermi-LAT HESS:2016 . The value $\gamma=3$ has been found in Tanaka:2008 for RX J1713.7-3946. At energies higher than 400 GeV, VERITAS and Fermi-LAT observe gamma-ray emission from Tychoâs SNR with power-law index $\sim 2.92$ Archambault:2017 .

To evaluate the parameter $\alpha$, general results for the particles’ spectral exponent $\gamma$, obtained in the framework of the diffusive shock acceleration theory extended to the case of anomalous transport with “Lévy flights” and “Lévy traps”, have been used (see [Lagutin A.A., 2017, to be published]). It was shown that for $\beta=0.8$ the spectral exponent $\gamma=2.8\div 2.9$ corresponds to nondiffusive transport with $\alpha\sim 1.7$.

The technique for determining the anomalous diffusivity was described in LagutinTymentcev:2004 . In this work, we used the following value $D_{0}\approx 1.5\cdot 10^{-3}$ pc1.7y*-0.8*.

As a result we have the cosmic ray spectrum

[TABLE]

A set of anomalous diffusion model parameters adopted in this paper is given in Table 1.

4 Results

In the framework of the model described above we have analysed the energy spectra and mass composition behaviour in the energy range up to the cut-off observed at $\sim 5\cdot 10^{19}$ eV. The energy spectra for all particles and elemental groups of nuclei are shown in Figure 1 compared with the data of different experiments. Along with good overall agreement with the representation of all the basic spectrum features, several key model predictions should be noted.

As already mentioned, the anomalous diffusion provides the knee in the energy spectrum as a pure propagation effect without any specific assumptions about maximum energy that a particle can achieve in a source. Another important feature of the model is that at the energy of the knee the spectral exponent, $\eta$, is equal to the injection spectra exponent $\gamma$. Thus, if anomalous diffusion is, indeed, the main mechanism of the knee feature, then the injection spectra exponent could be estimated from the analysis of the experimentally measured energy spectra around the knee.

The energies of knee retrieved for five elemental groups of nuclei (H, He, CNO, NeMgSi, Fe) are shown in Table 2. One can see that the knee energies are almost proportional to $Z$. We note that our estimations for the H spectrum are consistent with the results Bartoli:2015 where the knee-like feature in the (p+He) spectrum is observed at energies around $700\pm 230$ TeV.

Mean logarithmic mass and elemental fractions of CRs from $10^{15}$ to $10^{20}$ eV are shown in Figure 2. According to our model, the mass composition became heavier with energy up to $\sim 4\cdot 10^{17}$ eV. It is seen that at $10^{17.5}$ eV our data is consistent with LOFAR estimations Buitink:2016 .

The non-scattered particles of atomic charge $Z$ contribute to the total flux starting from the energy $3\cdot 10^{17}Z$ eV, at which their gyroradius became large enough to escape from the typical Galactic inhomogeneity and reach the observer. This results in a decreasing mean mass of CRs in the energy region ($4\cdot 10^{17}\div 2\cdot 10^{18}$ eV) matching the PAO data Settimo:2016 quite well. The lightest composition ($\langle\ln A\rangle\sim 1.6$) is observed at $\sim 2\cdot 10^{18}$ eV.

At ultra-high energies we have a gradually heavier composition with progressive cut-offs in each group of nuclei at $E_{0}=4\cdot 10^{18}Z$ eV. At $10^{19}$ eV our predictions give a heavier composition in comparison with the Auger data obtained from $X_{\max}$ distributions Settimo:2016 , but agree qualitatively with recent estimations Aab:2016 made on the basis of the advanced $X_{\max}$ and ground signal correlation analysis.

Finally it should be noted that our model contains an ankle in the light component similar to that observed by KASCADE-Grande Apel-a:2013 and more commonly the multiple ankle-like features in spectra of elemental groups of nuclei at $E>3\cdot 10^{17}Z$ eV. This should be considered as a crucial model prediction, which can be verified with the improved high-precision future measurements.

5 Conclusions

We have made a revision of the Galactic cosmic ray origin scenario in the framework of the anomalous diffusion model. It was shown that recent data on the cosmic ray energy spectrum and mass composition can be well reproduced by our model over the entire primary energy range with a reasonable set of key assumptions.

The knee energies were retrieved for different elemental groups of nuclei showing the rigidity-dependent behaviour similar to the standard scenario prediction. It is shown that, if the knee feature is caused by the anomalous diffusion of particles through the fractal-like ISM, the source injection spectra exponent could be estimated from the experimentally observed energy spectra being equal to the spectral exponent at the knee energy.

The following basic model predictions are expressed.

• •

The suppression of the all-particle spectrum at $E=5\cdot 10^{19}$ eV is due to the nuclei fluxes cut-offs caused by the rigidity-dependent energy limitation of the Galactic sources.

• •

The composition becomes heavier with energy from the knee to $\sim 10^{17.5}$ eV and reaches a maximum of mean logarithmic mass $\langle\ln A\rangle\sim 2.4$.

• •

In the energy region ($4\cdot 10^{17}\div 2\cdot 10^{18}$ eV) the mean logarithmic mass decreases reaching the minimum value of $\langle\ln A\rangle\sim 1.6$.

• •

The rapid weighting of the mass composition is observed at $E>2.5\cdot 10^{19}$ eV up to the pure iron composition at the cut-off.

• •

Multiple ankle-like spectrum features at energies above $3\cdot 10^{17}Z$ eV reflecting the contribution of the non-scattered CR component are expected.

Acknowledgements

This work was supported in part by the Ministry of Education and Science of The Russian Federation (state assignment for the fundamental and applied research performed at Altai State University).