Towards solving the mass-composition problem in ultra high energy cosmic rays

TL;DR

This paper uses shape-fitting of $X_{max}$ distributions to reconcile mass composition and energy spectra of ultra high energy cosmic rays across multiple experiments, proposing a dominant proton-helium model consistent with cosmology.

Contribution

It introduces a novel shape-fitting analysis of $X_{max}$ distributions that suggests a dominant p+He composition can explain observed spectra and composition data, considering the unreliability of proton-helium separation.

Findings

Proton-helium mixture can saturate the total flux in some models.

The p+He model explains the energy spectrum dip observed in recent data.

Potential for Auger-Prime to test the model's validity.

Abstract

Using the Auger mass-composition analysis of ultra high energy cosmic rays, based on the shape-fitting of $X_{max}$ distributions, we demonstrate that mass composition and energy spectra measured by Auger, Telescope Array and HiRes can be brought into good agreement. The shape-fitting analysis of $X_{max}$ distributions shows that the measured sum of proton and Helium fractions, for some hadronic-interaction models, can saturate the total flux. Such p+He model, with small admixture of other light nuclei, naturally follows from cosmology with recombination and reheating phases. The most radical assumption of the presented model is the assumed unreliability of the experimental separation of Helium and protons, which allows to consider He/p ratio as a free parameter. The results presented here show that the models with dominant p+He composition explain well the energy spectrum of the dip in the latest (2015 - 2017) data of Auger and Telescope Array, but have some tension at the highest energies with the expected Greisen-Zatsepin-Kuzmin cutoff. The Auger-Prime upgrade experiment has a great potential to reject or confirm this model.

Figures

Equations

$${\eta }_p(E)=\frac{J_p(E)}{J_p^{\mathrm{unm}}(E)},$$

$$\eta (E_{\mathrm{obs}})\pm \delta \eta =\frac{J_{\mathrm{obs}}(E_{\mathrm{obs}})}{K{E}_{\mathrm{obs}}^{-{\gamma }_g}}\times \left[1\pm \frac{1}{J_{\mathrm{obs}}(E_{\mathrm{obs}})}\frac{\partial J_{\mathrm{obs}}(E_{\mathrm{obs}})}{\partial E_{\mathrm{obs}}}\delta E_{\mathrm{obs}}\pm \frac{\delta J_{\mathrm{obs}}(E_{\mathrm{obs}})}{J_{\mathrm{obs}}(E_{\mathrm{obs}})}\right]$$

$$\begin{gathered} \frac{\delta (E^{3}J(E))}{E^{3}J(E)}=\frac{\delta J(E)}{J(E)} \\ +3\frac{\delta E}{E}, \end{gathered}$$

$$\frac{\delta \eta (E)}{\eta (E)}=\frac{\delta J(E)}{J(E)}+{\gamma }_g\frac{\delta E}{E}\text{leavevmode}.$$

Full text

1 Introduction

Mass composition still remains a controversial issue in Ultra High Energy Cosmic Rays (UHECR). The three biggest detectors, Pierre Auger (referred here as ’Auger’), Telescope Array (referred as ’TA’) and HiRes, have obtained contradictory results concerning mass composition of primary particles in the energy range (3 - 100) EeV (1 EeV = $1\times 10^{18}$ eV). At (1-3) EeV all three detectors agree with light composition, protons or protons and Helium, but in the range (3 - 100) EeV the Auger detector, the biggest one, founds a progressively heavier mass composition with increasing energy, while the other two detectors report the mass composition similar to that at lower energy.

At present there are two basic methods to study the mass composition of UHECR: direct measurements and indirect tests. The more reliable direct method is based on the observation of the fluorescent light produced by Extensive Air Showers (EAS) in the atmosphere. The indirect test is based on the signatures of mass composition in the primary energy spectrum and we start with it as more old and less constraining.

This approach works most efficiently for protons due to their interaction with the Cosmic Microwave Background (CMB). It results in two very specific spectral features: the Greisen-Zatsepin-Kuzmin (GZK) cutoff [2] and the pair-production dip. The former is a sharp cutoff at the end of the spectrum, around $E\sim 50$  EeV, due to the photo-pion production and the latter is a rather faint feature at $E\sim 1-30$ EeV first calculated in [3] and studied in detail in [4, 5, 6]. The dip is well confirmed in the spectra of all three detectors but its origin as the pair-production dip $p+\gamma_{cmb}\to p+e^{-}+e^{+}$ is now questioned by the Auger mass composition. Before 2011 the data published by HiRes [7] and Auger [8], and later confirmed by TA, showed high energy steepening in good agreement with the predicted GZK cutoff. Nevertheless, the newest data of Auger and TA, released in 2015 - 2017 seem to be in contradiction with this interpretation, see Fig.  1.

The propagation of UHE nuclei does not leave any clear signature of the mass composition in the energy spectra. The main channel of energy losses, that determines the spectrum of UHE nuclei, is the photo-disintegration process on the Extragalactic Background Light (EBL) and on CMB. This process naturally produces secondary lighter nuclei, mixing thus with the primary composition. As was first predicted by Gerasimova and Rozental [9] in 1961, i.e. before the discovery of CMB, nuclei photo-disintegration on EBL results in a suppression of the UHECR energy spectrum (GR steepening). In fact, as was realised later, see e.g. [10], a more sharp cutoff occurs at higher energies where the nucleus photo-disintegration time on CMB becomes equal to that on EBL. This cutoff arises at Lorenz-factor $\Gamma\sim(3-5)\times 10^{9}$ for all nuclei. The energy of the cutoff $E_{cut}\propto A\Gamma$ is different for primary nuclei with different $A$. This fact together with the unavoidable mixed composition, due to the production of secondary nuclei makes unclear any composition signature in the observed spectrum.

At present the best method to measure the mass composition is given by the observation of fluorescent light produced by the e-m component of EAS in the atmosphere. All three aforementioned detectors use this method. However, for better accuracy the fluorescent-light method needs additional information, which in the case of HiRes is given by the stereo observation of fluorescent light, and in the case of Auger (and recently of TA) this additional information is obtained from the data of on-ground detectors (water-Cherenkov detectors in Auger and scintillation detectors in TA).

The basic observable parameter related to mass composition is $X_{\max}(E)$, the atmospheric depth where the number $N(E)$ of particles in the cascade, with total energy $E$, reaches its maximum $X_{\max}$, is sensitive to the number of nucleons in the primary nucleus. Heavy nuclei interact higher in the atmosphere and have smaller fluctuations. In practice the actual quantity which allows to find the mass composition is the distribution $N(X_{\max})$ of the showers with total energy $E$.

In the case of large statistics the direct use of $N(X_{\max},E)$ gives the most reliable estimation of composition. In the case of limited statistics one may use the moments of this distribution, see e.g [11], namely the first moment which is the mean value $\langle X_{\max}\rangle$ and the second moment $\sigma(X_{\max})$ which is the variance or dispersion (RMS) of the distribution. As was demonstrated in [1] using only the first two moments for the analysis, may result in a false degeneracy: two different mass compositions may produce the same $\langle X_{\max}\rangle$ and $\sigma(X_{\max})$.

The shape-fitting analysis of $N(X_{max})$ recently performed by the Auger collaboration [1] gives very important results that, summarising, can be described as follows. The mass composition is assumed as a discrete sum of four elements: Iron (Fe), Nitrogen (N), Helium (He) and protons (p). For each element the $X_{\max}$ distribution is calculated by Monte Carlo (MC) simulations and the fraction of each element in the total flux is found from the comparison with observations. These fractions depend on the models of hadronic interaction included in MC simulations. A decisive result is given by the very small fraction of Iron at all energies, almost independently of the hadronic interaction model (see the upper panel of Fig. 4 in [1]). Besides, the analysis of [1] shows that the fraction of light elements (p+He) is quite large independently of the hadronic interaction model. It allows the conclusion that at least a large fraction of UHECR, if not the dominant one, is composed by light elements. The small fraction of Iron and large of p+He seem to be a common conclusion of the $N(X_{max})$ shape-fitting analysis of Auger and HiRes/TA data.

The argument above does not dismiss the question: Auger, HiRes and TA use the same fluorescent data to measure $X_{\max}$ and the same moments-based method for the analysis of mass composition. Why then their conclusions differ? The most convincing answer to this question is probably given in a recent paper by the TA collaboration [12]. The observation of fluorescence light can be performed in two ways: with a monocular observation, when only one telescope observes the fluorescent signal, or in the stereo mode when more than one telescope simultaneously observe the same shower. Fluorescence detection in monocular mode is less efficient to measure $X_{\max}$ in comparison with the stereo mode. HiRes and later TA used, apart from monocular, also stereo events with higher precision in the measurement of $X_{\max}$. It became possible because of the smaller (in comparison with Auger) spatial separation between telescopes. Auger, on the other side, to cover a larger area has a much larger separation among telescopes and collected mainly monocular fluorescent events. Instead, the Auger collaboration elaborated the hybrid technique, first proposed in [13], based on additional accompanying signal from, at least, one on-ground water-Cherenkov detector. Hybrid method allows to measure the core location and geometry of the shower, which improve the measurement precision for $X_{\max}$ and shower energy $E$. Auger collected now the largest number of hybrid events and we compare our predictions with the hybrid Auger data whenever this is possible.

At present TA is also using the hybrid technique with the help of 507 on-ground scintillation detectors [12]. With an accumulated statistics of 5 years data, TA reports [12] that the hybrid measurements of $X_{\max}$ agree with the results of Auger, if analysed with the EPOS-LHC hadronic interaction model [14]. On the other hand, using the QGSJetII-03 hadronic interaction model [15] the TA collaboration founds a mass composition compatible with only light nuclei.

Another important method to measure the mass composition is given by the observation of muons produced in the EAS.

The basic effect to distinguish a nucleus from a proton with the help of muons is related to the different energy per nucleon, $E/A$, at fixed total energy $E$. A low energy nucleon produces low energy charged pions which decay to muons before the parent pion undergoes new collisions with air-nucleus. Produced in the EAS, muons propagate rectilinearly with velocity $v\approx c$. As a result they can provide directional and timing informations, which can further reduce uncertainties in the fluorescent method. There are two well known muon quantities relevant for measuring mass composition: The total number of muons in the shower $N_{\mu}$ and the so-called Muon Production Depth (MPD) $X_{\max}^{\mu}$, which gives the atmospheric depth where the production rate of muons reaches its maximum [16, 17, 18].

The total number of muons $N_{\mu}$, called also the muon size of the shower, is especially important to determine mass composition at energies below the UHECR regime $<10^{17}$ eV, where the fluorescent emission is too faint to be detected. At these energies, $N_{\mu}$ gives the only way to determine the mass composition. The analysis of the muon component in the KASCADE-Grande experiment [19] allowed recently to find the Iron knee in the Galactic cosmic ray spectrum at energy $E\sim 1\times 10^{17}$ eV as predicted by the rigidity relation [20].

Among the three biggest UHECR arrays at present only the Auger experiment has several unique possibilities to measure the muon flux directly and use it to determine the mass composition. The on-ground water-Cherenkov detectors can measure muons in inclined directions, although with a high level of uncertainty due to the decoupling of the electron and muon signals. In AMIGA (Auger Muon and Infill Ground Array) there are muon detectors in the form of scintillation counters buried at a depth of 2.3m underground [21].

The new exciting method of muon detection in Auger experiment is given by the Auger-Prime [22] upgrade, recently funded, has been specifically designed to improve muon detection in the whole energy range of the experiment. Each water-Cherenkov tank will be equipped with scintillator layer on the top, sensitive only to e-m component of the shower, while water-Cherenkov detector is sensitive to both e-m and muon components. The combination of the two signals allows to reconstruct each of the fluxes separately. Recently, also the TA collaboration started important upgrades to increase the statistics at the highest energies enlarging the area covered by the surface detectors, an updated description of the status of the TA experimental set-up can be found in [23].

Another important channel of measuring the mass composition is connected with the correlation between muon signal and $X_{\max}$. We will discuss this method in some details in Section 4.

Above we discussed the problem of measuring the mass composition. Somewhat different question is whether it is possible to exclude, on the basis of observations, a pure proton composition at $E>3$ EeV. There are two challengers for such task: the quasi-isotropic gamma-radiation and neutrinos; both produced mainly by collisions of UHECR protons with background photons. The most stringent limit on the isotropic component of gamma-radiation in the range 50 MeV - 820 GeV is given by the Fermi-LAT experiment [24]. The strongest upper limit on the allowed UHE proton flux was obtained practically simultaneously in 2016 in three works [25, 26, 27]. The limit depends on the models for sources of UHECR, especially strongly on the injection power-law index $\gamma_{g}$ and cosmological evolution of sources. In [27] it was demonstrated that in a wide range of parameters the proton models which explain UHECR flux and spectrum, predict gamma-rays and neutrinos below the Fermi-LAT upper limit and IceCube flux.

In the present paper we use the latest Auger and TA observations, comparing them with the spectral features that arise due to propagation of UHECR, and their mass composition.

We argue that the spectral features may still be considered as an indication for a light mass composition, solving the problem of the alleged discrepancy between Auger and TA observations.

The paper is organised as follows: In Section 2 we reconsider the status of the dip model in light of the latest observations of the spectrum. In Section 3 we show how a mixture of Helium nuclei and protons provide a good description of the observed flux. In Section 4 we discuss the correlation of $X_{\max}$ with muon characteristics: $N_{\mu}$, the total number of muons in a shower, and $X_{\mu}^{\max}$, the Muon Production Depth; for this discussion we calculated the spectrum in the model p+He+CNO. The conclusion is given in Section 5.

2 Modification factor as indication of

proton-dominated composition

Propagating through CMB the proton energy spectrum acquires two characteristic features: the dip, due to the process $p+\gamma_{\rm cmb}\to p+e^{-}+e^{+}$ [3], and GZK cutoff [2] due to the reaction $p+\gamma_{\rm cmb}\to N+\pi$. These two features are quite different from the spectral features arising in the flux of UHE nuclei due to the interaction with CMB and EBL. This difference can be used to distinguish protons from nuclei and can be used as an additional (indirect) test of the mass composition, to be compared with other observations.

This test becomes particularly important in the light of the uncertainties in the direct measurements of the mass composition and in the hadronic interaction models at energies above the LHC (CERN) calibrations.

In this section we use the modification-factor method to identify protons in UHECR. Following the works [4, 5, 6], it can be proved that this method favours a proton-dominated mass composition in the observations of four experiments: AGASA, Yakutsk, HiRes and Auger, using the data before 2009, and the TA data as of 2011. Here we reconsider this analysis using the higher statistics data of Auger and TA as published in 2015-2017 [28, 29, 30, 31]. We will demonstrate that new data may agree at some reasonable conditions with the proton dip, but show noticeable differences with the GZK cutoff (see Fig. 1).

The proton modification factor used in the forthcoming analysis is defined as the ratio:

[TABLE]

where $J_{p}(E)$ is the total flux of protons, measured or calculated, taking into account all energy losses due to $p\gamma_{\rm cmb}$ collisions and adiabatic energy losses. In Eq.(1) we introduced also the unmodified proton spectrum $J_{p}^{\rm unm}(E)$ which is calculated taking into account only adiabatic energy losses, $J_{p}^{\rm unm}(E)=KE^{-\gamma_{g}}$.

Model-dependent phenomena enter both numerator and denominator in Eq.(1) and compensate or even cancel each other, while interaction with CMB photons does not enter the unmodified flux, i.e. the denominator, and appears only in the numerator of Eq.(1). Thus the modification factor presents, in an unsuppressed way, such features as dip and GZK cutoff directly connected with the propagation of UHE protons, while the model-dependent features are seen there in a suppressed form. Therefore the modification factor is an excellent instrument to search for the proton-dominated mass composition through the proton interaction features, dip and GZK cutoff, but it is not sensitive to the model-dependent features e.g. to the details of the acceleration models.

Modification factor depends non-trivially on the statistics of the events.

First, at very high observational statistics the agreement of the observed modification factor with the predicted one (both for protons) must be worse because at high resolution the ratio starts to distinguish better the model-dependent effects being less compensated in denominator and nominator of Eq.(1). In other words the higher statistics results in a worse agreement between denominator and numerator in Eq.(1) due to model-dependent effects. This phenomenon will be referred to as "high-statistics de-compensation".

The second effect is caused by the observation technique. There are two kinds of observational errors: flux errors $\delta J/J$ and energy errors $\delta E/E$. In all published presentations of the measured spectra the errors are given as the flux errors $\delta J/J$, while the energy errors $\delta E/E$ are just mentioned or shortly discussed. The systematic energy error was first presented in the Auger ICRC report 2013 [32, 33] as $\delta E/E=0.14$.

**In the present paper we consider two cases:

**

**(i) The total error is dominated by flux error $\delta J/J>>\delta E/E$, **

(ii) The energy error $\delta E/E$ is the dominant component $\delta E/E>\delta J/J$, or $\delta E/E>>\delta J/J$.

---

We start with the modification factor calculation $\eta(E_{\rm obs})\pm\delta\eta$ in the general case of two errors $\delta E_{\rm obs}$ and $\delta J_{\rm obs}$. It is easy to obtain the exact expression for $\eta(E_{\rm obs})$ in the case of natural condition $\delta E_{\rm obs}<<E_{\rm obs}$:

[TABLE]

In the case (i), $(\delta E/E)_{\rm obs}<<(\delta J/J)_{\rm obs}$, the second term in rhs of Eq.(2) may be neglected and the comparison of the proton modification factor with the data of four experiments released before 2011, namely AGASA, Yakutsk, HiRes, Auger, and also the 2011 data of TA, results in an excellent agreement with the observed spectra in approximately 100 energy bins [4, 5, 6]. It is a remarkable fact that this agreement is achieved with only one free parameter, the injection power-law index $\gamma_{g}\approx 2.6$, together with sources emissivity ${\cal L}_{0}$, which provides the total normalization of the flux using the same restriction $(i)$   $\delta E/E<\delta J/J$.

In Fig. 1 we compare the theoretical modification factor for protons with TA (left panel) Auger (right panel) data as of 2015 and 2017, using the same restriction $(i)\leavevmode\nobreak\ (\delta E/E)<\delta J/J$. This comparison shows a fairly good agreement with the dip in both datasets, the relatively small discrepancy with GZK cutoff for Auger and stronger discrepancy with GZK cutoff for TA. We postpone the discussion of GZK cutoff to a forthcoming publication.

As to more details concerning the dip, TA combined spectra of 2015 and 2017, show an excellent agreement with the theoretical modification factor, while Auger data show a good agreement with theoretical curve for 2015 hybrid data, but for the high-statistics 2017 combined data there is statistically significant difference with both the theoretical spectrum and the 2015 hybrid spectrum. Since the statistics of the combined Auger events is much higher than that for TA and for hybrid Auger events 2015, one may suspect the "statistical de-compensation"effect as the explanation for discrepancy, namely suppression of the compensation mechanism for model-dependent effects in numerator and denominator of modification factor given by Eq.  (1).

The reasonable agreement of the modification factor calculated above within the assumption $\delta E/E<<\delta J/J$ may be considered as experimental confirmation of this assumption, however there is an argument in favour of alternative relation (ii) $\delta E/E>\delta J/J$. From the Auger presentations of the combined energy spectra of 2014 - 2017 in the form $E^{3}J(E)$ with error-bar $\delta\left[E^{3}J(E)\right]$ :

[TABLE]

it is important to note that the total error $\delta(E^{3}J(E))$ is determined by the errors $\delta J$ and $\delta E$, which operate in two perpendicular directions. This remark is also valid for the modification factor $\eta(E)\propto E^{\gamma_{g}}J(E)$, which results in

[TABLE]

One may conclude that $\delta E/E$ term can be larger or much larger than $\delta J/J$ in the case of high statistics Auger combined events. Indeed, the first term $\delta J/J$ may be very small because of the tremendous Auger statistics (especially for the combined events) while the second term is large in particular for the systematic errors $\delta E/E=0.14$ [32]. Moreover, the difference between the hybrid and combined spectra in the Auger data (see below) could point toward a higher energy error. However, it is difficult to estimate it and we restrict our analysis to what discussed in [32]. Thus $\delta E/E\gg\delta J/J$, and it may over-close the whole dip in the vertical direction making it unobservable: the dip exists but cannot be seen because of too large energy errors. We will refer to this effect as "over-closing of gap".

This case is illustrated by Fig. 2 where the theoretical modification factor is plotted together with the Auger combined spectrum of 2017 and Auger hybrid spectrum of 2015. The two spectra show the largest difference in the recent Auger data. In the left panel they are presented with error bars $\delta J/J$, and in the right panel the energy errors $\delta E/E$ are added as half of systematic energy error of Auger [32, 33]. One may see that energy error-bars over-close the gap between the two Auger spectra and between each of Auger spectrum and theoretical modification factor (see captions to Fig. 2).

The conclusion above is valid for small energy errors $\delta E/E<\delta J/J$. In the case of larger energy errors $\delta E/E>\delta J/J$ the energy error $\delta E$ can exceed the gap between predicted and observed values, $\eta_{th}$ and $\eta_{obs}$, in particular at most disturbing energies $E\sim 4-6$ EeV, and contradiction between $\eta_{th}(E)$ and $\eta_{obs}(E)$ disappears.

In a more general discussion the good agreement of the proton theoretical modification factor with Auger and TA observations, at least at the dip energies $1-30$ EeV, is a strong indication of proton or proton-dominated composition. However, one cannot consider it as the the final proof. Indeed, on one side we observe the unique shape of the dip produced by $p+\gamma_{\rm cmb}\to p+e^{-}+e^{+}$ scattering. On the other side, in models with mixed nuclei composition, see for instance [34], one can obtain a theoretical spectrum with practically the same shape of the pair-production dip, but using more than 10 free parameters in the theoretical model. This result demonstrates that the very specific shape of the pair-production dip is not the unique one explaining the observed spectrum.

As discussed in the Introduction, there are a few other observations that can in principle challenge the proton-dominated mass composition: the Auger observations of mass composition using fluorescence and muon signals, the Fermi LAT data for the diffuse gamma-ray background and the IceCube data for astrophysical neutrinos. If all these experiments provide in future evidences against the proton-dominated composition one must conclude that the pair-production dip is an accidental coincidence. Nevertheless, until experimental data are not conclusive, one must consider the proton modification factor as an indication for a proton-dominated mass composition of UHECR.

3 $p+He$ model

As discussed in the Introduction there is some tension on the mass composition at $E>3$  EeV between the three biggest UHECR experiments Auger, TA and HiRes. In this Section we first summarise the basic physics of mass-composition measurements, then discuss the influence of cosmological environment and, finally, present calculations relevant to the $p+He$ model.

The measurement of mass composition is based on the $X_{\max}$ value, which is the depth of the atmosphere where the number of particles in the shower reaches its maximum. The value of $X_{\max}$ is a basic parameter to determine the mass composition of UHECR, with the best observable quantity for this determination given by the shape of the distribution $N(X_{\max},E)$ for showers with fixed total energy $E$.

As a matter of fact, until recently, instead of the distribution $N(X_{\max})$, the first two moments of this distribution were used: the mean value $\langle X_{\max}\rangle$ and its RMS $\sigma(X_{\max})$.

In two recent papers by the Auger collaboration, the mass-compositions obtained using moments-analysis ****[11]**** and shape-fitting $N(X_{max})$ analysis ****[1]**** are not identical. Realistically, they are not expected to be such, similarly to the already known fact that $\langle X_{\max}\rangle$ and $\sigma(X_{\max})$ give, if analysed separately, somewhat different results. The shape-fitting analysis is obviously the most fundamental and most sensitive method, since, for example, it involves the tiny parts of the wing distribution. Apart from it, the shape-fitting analysis demonstrated a degeneracy effect when two different mass compositions correspond to the same first two moments. For this reason, in the present paper, we choose the results obtained from the shape-fitting analysis of the Auger data ****[1]**** with the measured fractions of four nuclei species: Fe, N, He and p. These fractions, as determined from Auger measurements ****[1]****, reveal some uncertainties due to different hadronic interaction models namely QGSJet ****[15]****, EPOS-LHC ****[14]**** and Sybill ****[35]****.

The important result obtained in ****[1]**** is given by a very small fraction of Iron at all energies and for all interaction models, except EPOS-LHC at the two highest energy bins (see the upper panel in Fig. 4 of ****[1]****). The other important result, as mentioned in Introduction, and exposed in Fig. 3, is a large fraction of p+He, consistent with unity.

It is interesting to note that both effects have a natural cosmological explanation.

Among the heaviest nuclei, Iron is the most natural element to be produced in Supernova (SN) explosions and the absence of Iron in UHECR implies that other heavy elements must be absent too. Their suppression in the form $Fe/p\ll 1$ is very natural for extragalactic gas and extragalactic cosmic rays. Enhancement of p+He component has the same nature.

At the cosmological epoch of recombination, protons and Helium nuclei were the dominant components and heavy metals were almost completely absent. Meanwhile, the production of metals is compulsory. It is needed to provide cooling of ordinary stars during their evolution, including the preSN phase. The later stage of reionization in the universe, as detected by WMAP ****[36]**** and Planck ****[37]**** satellites occurs at redshift $z=11.0\pm 1.4$ and $z\lesssim 10.0$, respectively. This stage needs at least two early generations of stars with low metallicity, Pop III and Pop II stars. They inject into the extragalactic space a small amount of heavy metals. The main contribution to the Iron observed in the extragalactic space (and thus in extragalactic cosmic rays) is given by the present-time SN explosions. This scenario is confirmed by WMAP and Planck observations of the Universe reionization and by the observations of $Ly\alpha$ forest which indicate that extragalactic space had very low fraction of heavy elements at the level $Z\sim 10^{-3.5}Z_{\odot}$ at redshift $z\sim 5$, see eg. ****[38]****. Iron and other heavy metals are injected into extragalactic space mainly during a short interval $\Delta z$ at $z\sim 0$ mostly due to explosions of the last generation of SNe. This scenario is similar to the model of UHECR produced mostly nearby our Galaxy ****[39]****.

One may conclude that Hydrogen and Helium as the main products of primordial cosmological nucleosynthesis, and suppression of SN-produced Iron and other heavy metals in red-shifted gas, naturally result in a p+He dominated extragalactic gas and UHECR accelerated at red-shift z\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}afew.

In our calculations an additional simplifying assumption is used. Generically, we assume that all existing detectors do not distinguish reliably Helium from proton and one can consider p+He flux as one light component, assuming the fraction $He/p$ as a free parameter of the model. However, we will start with the sum p+He as it comes from Fig. 4 (strip 3 for He and strip 4 for p) of ****[1]****. Summing these two fractions, with errors summed in quadrature, we obtain p+He flux presented in Fig. 3. One can see that the sum of these two components saturates well the total flux, at least in the case of QGSJet and Sybil hadronic interaction models. This interesting fact confirms well our assumption that the light fraction (p+He) weakly depends on energy and with good accuracy saturates the total flux leaving small room for other components (e.g. $N$ which will be considered later as CNO component).

We are ready now to calculate the energy spectra for $p+He$ models and to compare them with spectra released by Auger and TA in 2015. We consider a power-law generation spectrum as $Q(E)=K_{i}E^{-\gamma_{g}}$ ($i=p,\leavevmode\nobreak\ He$) with the same generation index $\gamma_{g}$ for protons and Helium nuclei. We also assume that sources are distributed homogeneously and uniformly, so that the calculated spectrum is universal, i.e. not being affected by propagation models. Energy losses include pair-production, photo-pion production, and photo-dissociation for Helium on CMB and EBL.

Secondary protons from $He$ and $D$ photo-dissociation and also from neutron decays are included in calculations. For interaction with EBL photons the model ****[40]**** is used. In all these calculations we follow ****[34]****.

In Fig. 4 and Fig. 5 the computed spectra for $p+He$ models are presented for $\gamma_{g}=2.6$ and $\gamma_{g}=2.2$, respectively. A generic feature of $p+He$ spectra is the proton dominance at the highest energies because of the GR steepening for Helium at $E\lesssim 5\times 10^{18}$ eV due to photo-disintegration on EBL. Therefore the GZK cutoff in $p+He$ model becomes compulsory, unless the maximum acceleration energy $E_{\max}^{\rm acc}$ is below the GZK threshold $E_{\max}^{\rm GZK}\simeq 50$ EeV.

Consider first the case of generation index $\gamma_{g}=2.6$ shown in Fig. 4. This generation index corresponds to the canonical proton modification factor in the dip model (see section 2). Therefore if to take a small $He/p$ ratio at generation one should obtain the theoretical spectrum and theoretical modification factor in agreement with old (before 2015) observations. One may notice from Fig. 4 the similar agreement between theoretical and observational dips for Auger 2015 (hybrid data) and for TA 2015 (combined data). This agreement becomes worse at the highest energies.

In Fig. 4 we plot the comparison of the observed and calculated fluxes for Auger 2015 (left panel) and TA 2015 (right panel). It is remarkable that, at the dip energies, the TA spectrum can be described just rescaling by a factor $1.2$ the source emissivity needed to describe the Auger data. The behaviour of the flux at the highest energies is determined by the photo-pion production process. The maximum acceleration energy in Fig. 4 is taken at the level of $E_{max}=10^{21}$ eV. In other words the theoretical spectrum shape in Fig. 4 is exactly as predicted in the case of the GZK cutoff and, as follows from Fig. 4, it seems not well reproduced in both data sets. Auger shows an earlier cutoff at energies below the GZK cutoff energy ($\simeq 50$** EeV) while TA shows a flux suppression at energies slightly higher than this value.**

The fraction of Helium allowed at the sources depends on the assumptions for the injection power-law index. Assuming harder spectra it is possible to increase the fraction of Helium. In Fig. 5 we assume a maximum acceleration energy $E_{max}=8\times 10^{19}$ eV and a flatter injection spectrum with $\gamma_{g}=2.2$, that allows to increase the fraction of Helium nuclei in the generation spectrum up to $K_{He}/K_{p}=0.35$. This procedure improved but a little the agreement with observational data of Auger at the highest energies, while the good agreement with the dip remains practically as before. These changes are linked with the GR steepening of He spectrum due to photo-disintegration on the EBL radiation.

4 $X_{\max}$ and muons

As discussed in the Introduction, the observation of the showers muon component is a very efficient tool to measure the mass composition. The total number of muons $N_{\mu}$ in the shower and also the Muon Production Depth $X^{\mu}_{\max}$ are sensitive indicators for it. These two quantities are expected to have correlations with $X_{\max}$, measured using the fluorescent light, because all three quantities characterise the mass composition. Among the three biggest arrays, Auger, at present, is the most efficient one to measure the muon flux and the correlations of $N_{\mu}$ and $X^{\mu}_{\max}$ with $X_{\max}$.

To pursue such measurements the Auger collaboration has just started an overall upgrade of the experimental set-up (Auger-Prime) to instrument each water-Cherenkov detector with a $4m^{2}$ plastic scintillator located on the top ****[22]****, to disentangle the total signals produced by muons and electrons.

From the theoretical point of view, P. Younk and M. Risse ****[41]**** were the first to demonstrate that the statistical correlation of the shower-maximum depth $X_{\max}$ and the total number of muons $N_{\mu}$ in the shower, depends on the mixture of different nuclei species in the UHECR flux and thus can be used for the analysis of the mass composition. This is a transparent theoretical idea, because both quantities are good characteristics of the mass composition measured, however, in two independent experimental ways. In ****[41]****, the authors introduced the statistical correlation factor $r$, see Eq.(1) in ****[41]****, for the two quantities $X_{\max}$ and $N_{\mu}$ and demonstrated that $r$ is very sensitive to single nuclei composition and their mixtures. For example, a pure proton and a pure Iron composition gives $r_{p}=0.0$ and $r_{\rm Fe}=0.7$ respectively, while for equal ratios of both nuclei $r_{p+Fe}=-0.51$. These numbers are given for the case of an ’ideal detector’, for a realistic detector the difference is smaller. This example demonstrates the power of this method to shed additional light on the actual mass composition of UHECR.

Recently, restricting the analysis to hybrid events, the Auger collaboration presented the measurement of the correlation between the depth of shower maximum $X_{max}$, as observed by fluorescence telescopes, and the signal in the water-Cherenkov detectors ****[42]****. The surface array of these detectors is sensitive to muons, particularly in the case of inclined showers with zenith angles between $20$ and $90$ degrees, in which muons provide from $40\%$ to $90\%$ of the signal $S(1000)$ at a distance 1000 m from the shower core. Instead of $X_{\max}$ and $S(1000)$ in ****[42]**** as correlation quantities there were used $X_{\max}^{*}$ and $S_{38}^{*}$, which are the values of $X_{\max}$ and $S(1000)$ recalculated to a shower zenith angle of 38*∘*** and to a shower energy of 10 EeV. The statistical correlation factor measured and simulated in this way is denoted in [42] as $r_{G}(X_{\max}^{*},S_{38}^{*})$. The measured value obtained is $r_{G}^{\rm data}=-0.125\pm 0.024$ (negative), while the simulated values of the correlation factors for pure protons (p) and for a mixture of protons and Helium (0.8p+0.2He) were found approximately $r_{G}\approx 0$, for all three hadronic interaction models EPOS-LHC, QGSJetII-04 and Sibyll 2.1. This result seems to exclude both a pure proton model and a proton model with an admixture of 20% He, moreover, any p+He model, without an admixture of heavier nuclei, is disfavoured. The last conclusion is a strong argument against the pair-production dip, though this model may survive including He and more heavy nuclei as small admixtures.**

On the other hand one should realise that the present correlation measurement is undoubtedly a step forward to enter the muon physics in UHECR, since the correlation involves muon signals. Auger ****[43]**** and TA ****[44]**** face the problem here in the form of a muon excess by a factor 1.5 - 2.0 higher then expected at primary energies above 10 EeV. The analysis ****[45]**** shows that this excess is difficult to explain by any reasonable modification of the particle interactions and it maybe evidence for a light (proton) composition at these energies. The authors of ****[44]**** also interpret their measurements by a proton primary composition. Moreover, the muon excess questions the accuracy of the measurement of the muon signals based on inclined showers, and requires some caution towards the mass-composition restrictions obtained in ****[42]****. One may hope that in the near future the Auger-Prime measurements of $N_{\mu}$ and $X^{\mu}_{\max}$ will allow to find the correlation factors $r(X_{\max},N_{\mu})$ and $r(X_{\max},X^{\mu}_{\max})$ with the needed accuracy.

As discussed in the Introduction, a powerful method to determine mass composition with the help of muons is given by the measurement of the muon production depth of the shower $X_{\max}^{\mu}$ ****[16, 17, 18]****. Muons trace their parents, heavy nuclei or protons, being produced at different heights in the atmosphere and propagating rectilinearly with velocity $v\approx c$. The timing of muon signals according to detailed calculations ****[18]**** will improve the accuracy of $X_{\max}^{\mu}$ method as well as the accuracy of its correlation with $X_{max}$.

The problem with mass composition as it is found in ****[42]**** requires the presence of nuclei heavier than Helium. Here, assuming that the result of ****[42]**** will be confirmed in a more convincing way through the correlations $(X_{\max},X_{\max}^{\mu})$, and $(X_{\max},N_{\mu})$, we included in $p+He$ model also CNO nuclei with ratios at the source $He/p$=0.2 and $CNO/p$=0.1 (both allowed by Auger data ****[1]****). The spectra obtained are shown in figure 6, calculated following the computation scheme of ****[34]****, show a quite good agreement with experimental data of both Auger and TA. The latter, as before, are reproduced assuming a source emissivity multiplied by a factor 1.2 respect to the Auger case.

5 Conclusions

As far as mass composition is concerned there are three methods to analyse the fluorescence data: $\langle X_{\max}\rangle$,  $\sigma(X_{\max})$ and the shape-fitting analysis of $N(X_{max})$ distribution ****[1]****. As it is well known, the first two methods (moments of the $X_{\max}$ distribution) do not agree well between themselves and both disagree with the shape-fitting analysis (compare the mass composition obtained in ****[11]**** and in ****[1]****). In this paper we used the Auger shape-fitting analysis as the most reliable and free from false degeneracies, see ****[1]**** and Introduction.

In the shape-fitting analysis the mass composition is described in terms of four nuclei species: p, He, N (we consider it as CNO) and Fe (which can be considered as Iron group including the heavy metals). The results of this analysis are given as fractions of the fluxes of these four elements, which depend rather strongly on the hadronic interaction models used. The new and important result of this analysis is the very small fraction of Iron, compatible with zero, practically for all models of hadronic interactions. We argue that this result is natural for the standard cosmology with reionization of the universe.

Our first observation is that using QGSJet II-4 and Sybill 2.1 for the hadronic interaction model the sum of protons and Helium nuclei fractions saturates with good precision the total flux, while for EPOS-LHC it leaves more space for other elements especially at the lowest and highest energies. Thus a reasonable model could be a p+He dominated injection with a small admixture of CNO muclei.

Next we made the ad hoc assumption that at present all existing detectors cannot distinguish reliably Helium nuclei from protons and we calculated the spectra for Helium and protons considering them as a single component with the same injection power-law index $\gamma_{g}$ equal to 2.6 and 2.2 and taking the ratio He/p at the source to fit the spectra of Auger and TA. These ratios are 0.15 for $\gamma_{g}=2.6$ and 0.35 for $\gamma_{g}=2.2$. The calculated spectra are shown in Figs. 4 and 5.

The highest energy part of these spectra are always dominated by protons, because high-energy He nuclei are photo-disintegrated in collisions with EBL photons.

In the case of $\gamma_{g}=2.6$ the observed dip is mainly produced by protons: it is the canonical case of the dip model considered in section 2 in terms of the modification factor. The new element of this study is the inclusion of two errors: flux error $\delta J/J$ and energy error $\delta E/E$. At very large statistics of events the energy error dominates, and using $\delta J/J$ error e.g. for spectrum presented in the form $E^{3}J(E)$ is incorrect. It changes the status of the dip model in terms of the combined Auger spectrum (see Section 2).

In the case $\gamma_{g}=2.2$ the dip at EeV energies of the spectrum is produced by both Helium and proton components. However, the observed high energy cutoff in the Auger spectrum is located at energy below the predicted GZK cutoff. In any case the shape of the spectrum alone is not enough to accept the model (see Section 2).

Recently, the very interesting idea of the correlation of $X_{\max}$ and muon signals was proposed ****[41]**** and the Auger experiment ****[42]**** has developed it. Its results exclude a pure proton $p$ and $p+He$ mass composition. We argue in Section 4 that the accuracy of the found correlation is questioned by the anomalous muon excess observed in Auger ****[43]**** and TA ****[44]**** experiments. With near-future Auger-Prime measurements the accuracy of the measured muon signals $N_{\mu}$ and $X_{\max}^{\mu}$ will provide reliable values of the correlation functions $r(X_{\max},N_{\mu})$ and $r(X_{\max},X^{\mu}_{\max})$ together with the conclusions concerning mass composition. With these data the accuracy of the correlation analysis will undoubtedly reach the needed level providing important informations on mass composition.

In order to account for the correlation models of ****[42]****, considering them as preliminary, we excluded the pure two-component p+He models adding to them CNO nuclei. The CNO component seems also to appear in the analysis of ****[1]**** with the EPOS-LHC hadronic interaction model. As far as energy spectra are concerned, the model with ratios He/p=0.2 and CNO/p=0.1 gives better agreement with the observed spectrum than the two-component p+He model.

We conclude emphasising that the understanding reached so far on the mass composition of UHECR is still not conclusive. The observations of mass composition are still contradictory and cannot exclude a pure light composition, while the observations of spectra agree fairly well with such hypothesis. For these reasons the high energy muon program of Auger, especially the measurement of $X_{\max}^{\mu}(E)$, will be a crucial test of the models discussed in this paper.

We finish with the following note.

This paper is mainly focused on the impact of the Auger analysis ****[1]**** on the mass composition problem in UHECR. We avoided the discussion of some accompanying problems. In particular we calculated the spectra at $E\geq 1$ EeV to avoid the discussion about the transition from galactic to extragalactic cosmic rays. We also did not include the (possible) cosmological evolution of the sources, which can "artificially"improve the agreement of theoretical spectra with observations. We did not discuss the absence of GZK cutoff in all three biggest UHECR detectors: HiRes, Telescope Array and Auger. We hope to address these problems in a forthcoming publication.

The present paper is an updated version of the arXiv paper 1703.08671v1.

Acknowledgements

We thank Karl-Heinz Kampert, Antonio Villar, Piera Ghia, Michael Unger, Markus Risse and Alexey Yushkov for valuable remarks and critical comments. We are grateful to the anonymous Referee for many valuable remarks.

Список литературы

• [1] ** A. Aab et al. (Auger collaboration) Phys. Rev. D90, 122006 (2014). **
• [2] ** K. Greisen, Phys. Rev. Lett. 16 48 (1966).**

**G.T. Zatsepin and V.A. Kuzmin, JETP Letters 4 78, (1966). **

• [3] ** G.R. Blumenthal, Phys. Rev. D1 1596 (1970). **
• [4] ** V. Berezinsky, A.Gazizov, S. Grigorieva, Phys. Rev. D74 043005 (2006). **
• [5] ** V.Berezinsky, A.Gazizov, S. Grigorieva, Phys. Lett. B612 147 (2005). **
• [6] ** R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S. Grigorieva, B. Hnatyk, Astrop. Physics 27 76 (2007). **
• [7] ** T.Abu-Zayyad et al, Phys. Rev. Lett. 92, 151101 (2004). **
• [8] ** J. Abraham et al,Phys. Rev. Lett. 101, 061101 (2008). **
• [9] ** N.M. Gerasimova and I.L. Rozental, JETPh 41 488 (1961). **
• [10] ** R.Aloisio, V.Berezinsky, S.Grigorieva, Astropart. Phys. 41, 94, (2013). **
• [11] ** A. Aab et al. (Auger Collaboration), Phys. Rev. D90 122005 (2014). **
• [12] ** R.U.Abbasi et al (TA Collaboration), Astropart. Phys. 64 49 (2015). **
• [13] ** T. Abu-Zayyad et al. (HiRes-MIA collaboration), Astrophys.J. 557, 686 (2001). **
• [14] ** K. Werner, F.M.,Liu and T. Pirtoj, Phys. Rev. C74 044902 (2006). **
• [15] ** S. Ostapchenko, Phys. Rev. D74 014026 (2006). **
• [16] ** D. Garcia-Gamez et al. (Auger Collaboration), contribution ICRC 2011, arXiv:1107.4807. **
• [17] ** R. Conceicao, S. Andringa, L. Cazon and M. Pimenta (for the Auger collaboration), EPJ Web Conf. 52 03004 (2013), arXiv:1301.0507. **
• [18] ** L. Collica et al. (Auger Collaboration), Eur.Phys. J. Plus 131 (2016), arXiv:1609.02498. **
• [19] ** W.D. Apel et al (KASCADE-Grande Collaboration) Phys. Rev. Lett. 107 171104 (2011). **
• [20] ** V.S. Berezinsky, S.V. Bulanov, V.A.Dogiel, V.L.Ginzburg, V.S. Ptuskin, Astrophysics of Cosmic Rays, North Holland (1990). **
• [21] ** B. Wundheiler et al. (Auger Collaboration), contribution ICRC 2015, PoS(ICRC2015) 324 (2015), arXiv:1509.03732. **
• [22] ** A. Aab et al. (Auger Collaboration), The Pierre Auger Observatory Upgrade - Preliminary Design Report, arXiv:1604.03637. **
• [23] ** R.U. Abbasi et. al. (TA Collaboration), Astrophys.J. 858 (2018) no.2, 76. **
• [24] ** A.A. Abdo et al (Fermi-LAT Collaboration), Phys. Rev. Lett. 104 101101 (2015). **
• [25] ** R.Y. Liu, A.M. Taylor, X.Y. Wang, F. Aharonian, Phys. Rev.  D94 043008 (2016). **
• [26] ** E.Gavish, D.Eichler, Astrophys.J. 822 56 (2016). **
• [27] ** V. Berezinsky, A. Gazizov, O. Kalashev Astropart. Phys. 84 52 (2016). **
• [28] ** I. Valino et al. (Auger Collaboration), contribution ICRC 2015, PoS (ICRC2015) 271 (2015), arXiv:1509.03732. **
• [29] ** D. Ivanov et al. (TA Collaboration), contribution ICRC 2015, PoS (ICRC2015) 349 (2015). **
• [30] ** F. Fenu et al. (Auger collaboration), contribution ICRC 2017, PoS (ICRC2017) 486 (2017), arXiv:1708.06592 **
• [31] ** J. Matthews et al. (TA collaboration), contribution ICRC 2017, PoS (ICRC2017) 1096 (2017). **
• [32] ** V. Verzi et al. (Auger collaboration), contribution ICRC 2013, arXiv:1307.5059. **
• [33] ** V. Verzi, rapporteur report ICRC 2015, PoS ICRC2015 015 (2016). **
• [34] ** R. Aloisio, V. Berezinsky, P. Blasi, JCAP 10 020 (2014). **
• [35] ** E.J. Ahn et al, Phys.ReV. D80 094003 (2009). **
• [36] ** G. Hinshaw et al. (WMAP Collaboration) Astroph. J. Suppl. Ser. 180 225 (2009). **
• [37] ** R. Adam et al. (Planck Collaboration), Astron. Astrophys. 596 A108 (2016). **
• [38] ** A. Songaila, Ap.J. 561 L153 (2001). **
• [39] ** R.Y. Liu, A.M. Taylor, X.Y. Wang, F.A. Aharonian Phys.Rev. D94 (2016) no.4, 043008, arXiv:1603.03223. **
• [40] ** F.W. Stecker, M.M. Malkan and S. Scully, Astrophys. J. 648 774 (2006). **
• [41] ** P. Younk and M. Risse, Astrop. Phys. 35 807 (2012). **
• [42] ** A. Aab et al. (Auger Collaboration) Phys. Lett. B708 288 (2016). **
• [43] ** A. Aab et al for Auger collaboration, Phys. Rev. Lett. 117 192001 (2016). **
• [44] ** R. Takcisin et al. JPS Conf. Proc. 19 011045 (2018). **
• [45] ** S. Ostapchenko, XXV ECRS 2016 Proceedings, arxiv:1612.09461. **