text:
# About EAS size spectra and primary energy spectra in the knee region

## 1 Introduction

High statistical accuracy in modern EAS experiments in the knee region encouraged the investigation of the fine structure of EAS size spectra. Although the origin of the knee is still a matter of debate, recently the series of publications  appeared, where the sharpness and the spectral structure in the knee region were interpreted by the contribution of heavy nuclei from a single local supernova. Along with this, the absolute differential EAS size spectra measured at different atmosphere depths and different zenith angles are not explained yet from the point of view of a single $`AA_{Air}`$ interaction model and a single model of primary energy spectra and elemental composition. Such an attempt has been made in  based on the QGS interaction model and rigidity-dependent steepening primary energy spectra  for a description of vertical MSU, AKENO and Tien-Shan EAS size spectra. In the present work we worked out a formalism of the inverse problem solution - reconstruction of the primary energy spectrum and elemental composition based on the known EAS size spectra of KASCADE , AKENO , EAS-TOP  and ANI  measured at different zenith angles. The calculations were done in the frames of QGSJET  and SIBYLL  interaction models. As a primary spectrum we have tested the modified rigidity-dependent steepening primary energy spectra and the hypothesis of the additional component in the knee region . In this case the type of nucleus of the additional component was considering as unknown and determining by the best fit of the fine structure of EAS size spectra in the knee region.

## 2 EAS inverse problem

In general, the energy spectra ($`\mathrm{}_A/E_0`$) of primary nuclei ($`A`$) and detectable EAS size spectra ($`I/N_e^{}`$) are related by the integral equation -

$$\frac{I(E_e,\overline{\theta },t)}{N_e^{}}=\underset{A}{}_{E_{min}}^{\mathrm{}}\frac{\mathrm{}_A}{E_0}W_\theta (E_0,A,N_e^{},\overline{\theta },t)𝑑E_0$$
(1)
where $`E_0,A,\overline{\theta }`$ are energy, nucleon number (1-59) and average zenith angle of primary nuclei, $`E_e`$ is an energy threshold of detected EAS electrons, $`N_e^{}(E>E_e)`$ is the estimation value of EAS size obtained by the electron lateral distribution function. Here, by the EAS size ($`N_e(E>0)`$) we mean the total number of EAS electrons at given observation level $`(t)`$. The kern ($`W_\theta `$) of integral equation (1) is determined as

$$W_\theta \frac{1}{\mathrm{\Delta }_\theta }_{\theta _1}^{\theta _2}_0^{\mathrm{}}\frac{\mathrm{\Omega }(E_0,A,\theta ,t)}{N_e}P_\theta \frac{\mathrm{\Psi }(N_e)}{N_e^{}}\mathrm{sin}\theta d\theta dN_e$$

where $`\mathrm{\Omega }/N_e`$ is an EAS size spectrum at the observation level $`(t)`$ for given $`E_0,A,\theta `$ parameters of a primary nucleus and depends on $`AA_{Air}`$ interaction model; $`\mathrm{\Delta }_\theta =\mathrm{cos}\theta _1\mathrm{cos}\theta _2`$;

$$P_\theta P(Ne,E_0,A,\theta )=\frac{1}{XY}D(N_e,E_0,A,\theta ,x,y)𝑑x𝑑y$$

is a probability to detect an EAS by scintillation array at EAS core coordinates $`|x|<X/2`$, $`|y|<Y/2`$ and to obtain the estimations of EAS parameters ($`N_e^{}`$, $`s`$ \- shower age, $`x^{},y^{}`$ \- shower core location) with given accuracies; $`\mathrm{\Psi }/N_e^{}`$ is a distribution of $`N_e^{}(N_e,s,x,y)`$ for given EAS size ($`N_e`$).
One may achieve significant simplification of equation (1) providing the following conditions during experiments:
a) selection of EAS cores in a range where $`P_\theta 1`$,
b) the log-Gaussian form of the measuring error ($`\mathrm{\Psi }/N_e^{}`$) with an average value $`\mathrm{ln}(N_e\delta )`$ and a RMSD $`\sigma _N`$, where $`\delta `$ involves all transfer factors (an energy threshold of detected EAS electrons, $`\gamma `$ and $`\mu `$ contributions) and slightly depends on $`E_0`$ and $`A`$,
c) transformation (standardization) of the measured EAS size spectra to the EAS size spectra at observation level

$$\frac{I(0,\overline{\theta },t)}{N_e}\eta \frac{I(E_e,\overline{\theta },t)}{N_e^{}},$$

where $`\eta =\delta ^{(\gamma _e1)}\mathrm{exp}\{(\gamma _e1)^2\sigma _N^2/2\}`$ and $`\gamma _e`$ is the EAS size power index,
d) consideration of either all-particle primary energy spectrum $`\mathrm{}_\mathrm{\Sigma }/E_0`$ with effective nucleus $`A_{eff}(E_0)`$ or energy spectra of primary nuclei ($`\mathrm{}_{A_\xi }/E_0`$) gathered in a limited number of groups ($`\xi =1,\mathrm{}\xi _{\mathrm{max}}`$) as unknown functions.
(a-d) conditions make EAS data of different experiments more comparable and equation (1) converts to the form

$$\frac{I(0,\overline{\theta },t)}{N_e}=\eta _{E_{min}}^{\mathrm{}}\frac{\mathrm{}_\mathrm{\Sigma }}{E_0}\frac{\mathrm{\Omega }(E_0,A_{eff}(E_0),\overline{\theta },t)}{N_e}𝑑E_0$$
(2)
or

$$\frac{I(0,\overline{\theta },t)}{N_e}=\eta \underset{\xi =1}{\overset{\xi _{\mathrm{max}}}{}}_{E_{min}}^{\mathrm{}}\frac{\mathrm{}_{A_\xi }}{E_0}\frac{\mathrm{\Omega }(E_0,\overline{A}_\xi ,\overline{\theta },t)}{N_e}𝑑E_0$$
(3)
However, even in this form the determination of primary energy spectra by measured EAS size spectra and solution of integral equations (2,3) in general is unsolvable problem. At the same time, using the a priori information about energy spectra of primary nuclei ($`\mathrm{}_{A_\xi }/E_0`$) and the EAS size spectra $`I/N_e^{}f_{i,j}(N_{e,i}^{},\overline{\theta }_j,t)`$ measured in $`i=1,\mathrm{}m`$ size intervals and $`j=1,\mathrm{}n`$ zenith angular intervals, one may transform the inverse problem into $`\chi ^2`$-minimization problem

$$\mathrm{min}\{\chi ^2\}\mathrm{min}\left\{\underset{i}{\overset{m}{}}\underset{j}{\overset{n}{}}\frac{(f_{i,j}F_{i,j})^2}{\sigma _f^2+\sigma _F^2}\right\}$$
(4)
with unknown (free) spectral parameters. Here $`F_{i,j}F(N_{e,i}^{},\overline{\theta }_j,t)`$ are the expected EAS size spectra determined at the right hands of equations (1-3) and $`\sigma _f,\sigma _F`$ are the uncertainties (RMSD) of measured ($`f_{i,j}`$) and expected ($`F_{i,j}`$) shower size spectra.
One may also unify the data of different experiments applying minimization $`\chi _U^2`$ with re-normalized EAS size spectra

$$\mathrm{min}\{\chi _U^2\}\mathrm{min}\left\{\chi ^2(\frac{f_{i,j,k}}{_i_jf_{i,j,k}},\frac{F_{i,j,k}}{_i_jF_{i,j,k}})\right\}$$
(5)
where index $`k=1,\mathrm{}l`$ determines the observation levels ($`t`$) of experiments. Expression (5) offers an advantage for experiments where the values of methodical shift ($`\delta `$) and measuring error ($`\sigma _N`$) are unknown or known with insufficient accuracy.
The energy spectra of primary nuclei are preferable to determine (a priori) in the following generalized form

$$\frac{\mathrm{}_A}{E_0}\beta \mathrm{\Phi }_AE_0^{\gamma _1(A)}\left(1+\left(\frac{E_0}{E_{knee}(A)}\right)^ϵ\right)^{(\gamma _1(A)\gamma _2)/ϵ}$$
(6)
Unknown (free) spectral parameters in approximation (6) are $`\beta `$, $`E_{knee}(A)`$ (so called ”knee” of energy spectrum of $`A`$ nucleus), $`\gamma _1`$ and $`\gamma _2`$ (spectral asymptotic slopes before and after knee), $`ϵ`$ (sharpness parameter of knee, $`1ϵ10`$). The values of $`\mathrm{\Phi }_A`$ and $`\gamma _1(A)`$ parameters are known from approximations of balloon and satellite data  at $`A1,4,\mathrm{}59`$ and $`E_0110^3`$ TeV. Parameter $`\beta 1`$ determines the normalization of spectra (6) in $`10^210^5`$ TeV energy range.
Thus, minimizing $`\chi ^2`$\- functions (4,5) on the basis of measured values of $`I(\overline{\theta }_{i,k})/N_{ej,k}`$ and corresponding expected EAS size spectra (2,3) at given $`m`$ zenith angular intervals, $`n`$ EAS size intervals and $`l`$ experiments one may evaluate the parameters of the primary spectrum (6). Evidently, the accuracies of solutions for spectral parameters strongly depend on the number of measured intervals ($`mnl`$), statistical errors and correctness of $`\mathrm{\Omega }(E_0,A,\theta ,t)/N_e`$ determination in the framework of a given interaction model.

## 3 Results

Here, the parametric solutions of the EAS inverse problem are obtained on the basis of KASCADE  ($`t=`$1020 g/cm2), AKENO  (910 g/cm2), EAS-TOP  (810 g/cm2) and ANI  (700 g/cm2) published EAS size spectra. These experiments were carried out at different observation levels and were chosen for two reasons: satisfaction of (a-c) conditions from the section (2) and a high statistical accuracy of presented data (especially KASCADE experiment). Unfortunately, during the standardization of EAS size spectra (condition (c)) the value of $`\eta `$ parameter is not always known with proper accuracy, which is the main reason of discrepancy in the results of different experiments. In our calculations we included $`\eta `$ in the list of unknown spectral parameters and determined by the minimization of functional (4). The problem does not exist if the minimization of re-normalized EAS size spectra is unified, since the linear parameters ($`\eta _k\beta `$) are canceled out from the functional (5).
The differential EAS size spectra $`\mathrm{\Omega }(E_0,A,\theta ,t)/N_e`$ for given $`E_00.032,0.1,\mathrm{},100`$ PeV, $`A1,4,12,16,28,56`$, $`t0.5,0.6,\mathrm{},1`$ Kg/cm2, $`\mathrm{cos}\theta 0.8,0.9,1`$ were calculated using CORSIKA562(NKG) EAS simulation code  at QGSJET  and SIBYLL  interaction models. Intermediate values are calculated using 4-dimensional log-linear interpolations. The estimations of errors of the expected EAS size spectra $`\mathrm{\Omega }/N_e`$ at fixed $`E_0,A,\theta ,t`$ parameters did not exceed $`35\%`$.
The basic results of minimizations (4,5) at a given number ($`\nu `$) of unknown spectral parameters and the values of $`\chi ^2/q`$ (or $`\chi _U^2/q_u`$ for unified data), are presented in Tables 1-4, where $`q=mn\nu 1`$ and $`q_u=_k(mn)_k\nu 1`$ are corresponding degrees of freedom. The upper (lower) rows of each experiment in Tables 1 and 4 correspond to parameters obtained by the QGSJET (SIBYLL) interaction model.

### 3.1 Test of rigidity-dependent energy spectra

Table 1 contains the approximation values of spectral parameters at rigidity-dependent approach

$$E_{knee}(A)=RZ$$
(7)
where $`Z`$ is a charge of $`A`$ nucleus and R is a parameter of magnetic rigidity (or a critical (cutoff) energy in more modern model ). The results were obtained by minimizations (4,5) applying (6,7). The magnitudes of $`\gamma _1(A)`$ for all nuclei were taken from . The number of unknown (free) parameters is equal to $`\nu =4`$ and corresponding solutions at two interaction models are presented in Table 1.

The stability of solutions (or the steep of $`\chi ^2`$ global minimum) of minimizations (4,5) is seen from obtained errors. So, it is seen that the $`\chi ^2`$ minimum for AKENO data does not depend on the $`\epsilon `$ sharpness parameter at SIBYLL and partly at QGSJET interaction models.
The obtained slopes of primary energy spectra after knee agree with the same calculations  performed by QGS model and exceed well known expected values ($`33.1`$) in the $`10^{17}`$ eV energy range . Such a steep of primary spectra after the knee is a result of $`(\gamma _2,R)`$ correlations in (4-7). In case of $`R2000`$ TeV the $`E_k(Fe)5.210^4`$ TeV and the primary energy spectra in the large interval ($`E_k(H)E_k(Fe)`$) at fixed $`\gamma _1(A_\xi )`$ can be conformed with corresponding EAS size spectra provided abnormally steep slope ($`\gamma _23.43.5`$) after the knee.
The results of expected EAS size spectra in comparison with corresponding experimental data are shown in Fig. 1 (the thin solid lines by QGSJET model, the thin dashed lines by SIBYLL model). Despite the satisfactory agreement ($`\chi ^21`$) of EAS size spectra with predictions of rigidity-dependent steepening spectra (6,7) and QGSJET interaction model the form (fine structure) of the measured EAS size spectra in the knee region differs from the form of corresponding expected spectra. It is worth mentioning that the difference is formally small and does not exceed several percents in precise KASCADE data.
As a next step, we attempted to test the rigidity-dependent approximation (7) directly. The knee-energies $`E_k(A_\xi )`$ were chosen as unknown (free) parameters in the approximations of EAS data (4-7), where $`\xi =1,\mathrm{}\xi _{\mathrm{max}}`$. However, the stable solutions for free parameters were obtained only at unified minimization (5), fixed values $`\gamma _1(A_\xi ),\gamma _2,\epsilon `$ and number of nuclear groups $`\xi _{\mathrm{max}}5`$. In Table 2 the values of spectral parameters $`E_{knee}(A_\xi )`$ obtained by minimization of $`\chi _U^2`$ (5) for 5 groups of primary nuclei ($`H`$), ($`He,Li`$), ($`BeNa`$), ($`MgCl`$), ($`ArNi`$) at given values of $`\gamma _2=3.1`$ , $`ϵ=4.0`$ and $`\nu =5`$ are presented. It is seen, that approach (7) is performed only for nuclei with $`A>1`$ and $`R400`$ TV or (7) is valid for all nuclei at $`R400`$ TV but there is an additional proton flux with $`E_{knee}^{(p)}3000`$ TeV which shifts the knee value of the total proton energy spectrum.

The following testing of the rigidity-dependent approach (7) is based on the investigation of the all-particle energy spectrum. Toward this end the fits of primary energy spectra $`d\mathrm{}_A/dE_0`$ for different nuclei ($`A=1,\mathrm{}59`$) known from  were extrapolated up to $`E_A=10^5`$ TeV energies taking into account (6,7) at $`R600`$ TeV . The obtained expected all-particle energy spectrum $`dI_\mathrm{\Sigma }/dE_0=_Ad\mathrm{}_A/dE_0`$ was approximated by expression similar to (6) at five free parameters ($`\mathrm{\Phi }_\mathrm{\Sigma },\gamma _1,\gamma _2,E_k,\epsilon `$). The average values of primary nuclei ($`\mathrm{exp}(\overline{\mathrm{ln}A})`$) obtained from extrapolations of spectra  were approximated by step function

$$A_{eff}(E_0)=a+b\mathrm{ln}\left(\frac{E_0}{E_k}\right)$$
(8)
where $`b=b_1`$ at $`E_0<E_k`$ and $`b=b_2`$ at $`E_0>E_k`$. The results of these approximations are presented in a last row of Table 3 . The value of sharpness parameter was equal to $`\epsilon =1\pm 0.1`$.
The first and second rows of Table 3 content the parameters of the all-particle energy spectra ($`\mathrm{}_\mathrm{\Sigma }/E_0`$), which were obtained by minimization $`\chi _U^2`$ (expressions 2,5-7) of unified EAS size data at $`\nu =6`$ and $`\epsilon =1`$. The approximation (8) has been used for $`A_{eff}(E_0)`$ which is at the right hand of expression (2). It is necessary to note that the solutions for $`a,b_1,b_2`$ parameters can be obtained only by re-normalized EAS size spectra (5) because of the strong correlation between a linear parameter $`\eta `$ and an effective nucleus $`A_{eff}`$.
It is seen from Table 3 that the model of rigidity-dependent energy spectra predicts the increase ($`b_1>0`$, third row of Table 3) of the mean nucleus with energy, whereas the presented analysis of EAS data points out a decrease of $`A_{eff}`$ with energy ($`b_1<0`$, first two rows of Table 3) in the energy range of $`E_0<E_k`$. It is obvious, that despite $`\overline{A}`$ could not be exactly equal to $`A_{eff}`$, at least their dependence on energy must be the same.
The results of recent precise experiments DICE  and CASA-BLANKA  also point out to the decrease of $`\overline{\mathrm{ln}A}`$ with energy at $`E_0<E_k`$. This dependence of $`A_{eff}`$ and $`\overline{\mathrm{ln}A}`$ on energy might be explained by the contribution of an additional light component in a primary nuclei flux.

From the above analyses follows that rigidity-dependent steepening energy spectra in combination with QGSJET or SIBYLL interaction model can not explain the obtained results of the fine structure of EAS size spectra  in the knee region (Table 1 and Fig. 1), the large values of knee $`E_k(H)`$ for Hydrogen component (Table 2) and dependence of the effective nucleus $`A_{eff}(E_0)`$ on primary energy before the knee (Table 3). In this connection we have carried on the search of more adequate model of primary energy spectra and elemental composition.

### 3.2 Test of additional component

Based on predictions  the primary energy spectra in approximation (6) have been added by a new (polar cap ) component $`\mathrm{}_{Add}/E_0`$ with power energy spectrum

$$\frac{\mathrm{}_{Add}}{E_0}=\mathrm{\Phi }^{(p)}\left(E_k^{(p)}\right)^{\gamma _1^{(p)}}\left(\frac{E_0}{E_k^{(p)}}\right)^{\gamma ^{(p)}}$$
(9)
where $`\gamma ^{(p)}=\gamma _1^{(p)}`$ at $`E_0<E_k^{(p)}`$ and $`\gamma ^{(p)}=\gamma _2`$ at $`E_0>E_k^{(p)}`$.
New spectral parameters $`\mathrm{\Phi }^{(p)},\gamma ^{(p)},E_k^{(p)}`$ and nucleon number $`A^{(p)}`$ of the additional component are considered as unknown and determined together with parameters of rigidity-dependent energy spectra (6,7) by minimization of $`\chi ^2`$ and $`\chi _U^2`$ (4,5) at $`\nu =7`$. The results of expected EAS size spectra for each experiment (KASCADE, AKENO, EAS-TOP and ANI) taking into account contribution of additional component (9) are shown in Fig. 1 (the thick solid lines by QGSJET model and the thick dashed lines by SIBYLL model). It is seen that the additional component with high accuracy (2-5% for KASCADE and ANI data) describes the fine structure of EAS size spectra in the knee region. The values of slopes of additional component before the knee turned out to be $`\gamma _1^{(p)}1.5\pm _{0.2}^{0.6}`$. The nucleon number ($`A^{(p)}`$) of this component with high reliability did not exceed of $`A^{(p)}=1`$ for most of experiments especially at QGSJET interaction model (except from AKENO ($`A^{(p)}56`$)). The unified analyses of all experiments at QGSJET and SIBYLL interaction models also gave a proton or neutron ($`A^{(p)}=1`$) composition of the additional component. The values of other spectral parameters at $`\gamma _1^{(p)}=1.5`$ and $`A^{(p)}=1`$ are included in Table 4.
The obtained result disagrees with  where the alike component consists of several heavy nuclei. However, our result is based on the high accuracy of the coincidence ($`23\%`$ for KASCADE data) of expected and measured EAS size spectra by both $`\chi ^2`$ criterion and the overlapping of fine structure of spectra.
The comparison of parameters $`\gamma _2`$ from Tables 1,4 shows that spectral slopes after the knee from Table 4 roughly overlap with the expected slope ($`33.1`$) well known from $`N_e10^7`$ EAS data. This can be explained the fact that the spectral brake of summary proton component ($`E_k^{(p)}210^3`$ TeV) is closer to the iron component ($`E_k(Fe)10^4`$ TeV) at $`R400`$ TV.
The coexistence of rigidity-dependent primary energy spectra and additional proton flux with spectral parameters $`E_k^{(p)}2000`$ TeV, $`\gamma _1^{(p)}1.5`$, $`\gamma _2^{(p)}=\gamma _23.1`$ explains also the shifted value of the knee $`E_k(H)`$ in Table 2 and the decrease of $`A_{eff}(E_0)`$ at $`E_0<E_k`$ (see section 3.1).
Although the additional component does not contribute to shower sizes below or above the knee significantly , it is essential for the reproducing of the sharp knee of EAS size spectra. This idea is taken from  and is confirmed here.

The final all-particle energy spectrum ($`I/E_0`$) obtained by unified EAS size data at QGSJET interaction model

$$\frac{I}{E_0}=\beta \left(\underset{A}{}\frac{\mathrm{}_A}{E_0}+\frac{\mathrm{}_{Add}}{E_0}\right)$$

and corresponding energy spectra ($`\mathrm{}_A/E_0`$) of 6 nuclear groups with additional component ($`\mathrm{}_{Add}/E_0`$) at normalization $`\beta =1`$  are presented in Fig. 2. The solid (dashed) line is the all-particle energy spectrum obtained by unified (only KASCADE) EAS size spectra. Dotted lines are the energy spectra of different primary components obtained by unified EAS data. Symbols in Fig. 2 are the data from DICE , CASA-BLANKA  arrays and reviews .
The numerical values of all-particle energy spectrum (solid line in Fig. 2), corresponding error, spectrum of additional proton component and average nucleus versus primary energy are presented in Table 5.

## Conclusion

High statistical accuracy of experiments KASCADE, EAS-TOP, AKENO and ANI allowed to obtain approximations of primary energy spectra and elemental composition with accuracy $`15\%`$ in the knee region. Along with this, KASCADE and ANI EAS size spectra at $`\theta <37^0`$ are described with the accuracy of $`25\%`$ in a whole measurement interval.
Obtained results show the evidence of QGSJET interaction model at least in $`10^510^7`$ TeV energy range, rigidity-dependent steepening primary energy spectra at $`R200400`$ TV and existence of the additional proton (or neutron) component with spectral power index $`\gamma _1^{(p)}1.5\pm _{0.2}^{0.6}`$ before the knee $`E_k^{(p)}2030`$ TeV. The contribution of the additional proton (neutron) component in all-particle energy spectrum turned out to be $`20\pm 5\%`$ at primary energy $`E_0=E_k^{(p)}`$.

## Acknowledgements

We thank Peter Biermann for extensive discussions and Heinigerd Rebel, Johannes Knapp and Dieter Heck for providing the CORSIKA code.