Universality of Electron Distributions in Extensive Air Showers

Cosmic rays of the highest energies, ${E}_{0}\gt {10}^{17}\,\mathrm{eV}$, can be detected and studied only by registering the effect they produce in the atmosphere, i.e., extensive air showers, cascades of particles (mainly electrons of both signs), and photons reaching the ground. Thanks to the fluorescence technique, it is also possible to gain information about the state of a shower developing in the atmosphere at altitudes higher than the ground. Shower particles cause the atmosphere to emit fluorescence and Cherenkov light, the intensity and spectrum of which depend on the electron state at the time of emission. Thus, when measured many different times while the shower traverses the atmosphere, the light characteristics enable the derivation of the total energy deposit of the shower particles changing with atmospheric depth. This provides the best estimate of the shower primary energy, i.e., the energy of the cosmic-ray particle E₀. It also allows one to calibrate showers, which are registered by ground detectors all day, much more frequently and not only at dark nights, as is necessary for light detector measurements. This was done for the first time at the Pierre Auger Observatory (Aab et al. 2015).

However, an exact determination of the shower atmospheric profile, from which E₀ and X_max, the depth of the shower maximum, can be derived demands detailed knowledge of the state of the shower of the most abundant particles—electrons—at any stage of its development. Our present study shows that it is possible to obtain this profile, in contrast to the common belief that extensive air showers are a highly fluctuating phenomenon. Indeed, the atmospheric depths of cosmic particle interactions, the energies transferred to the secondary particles, and the consequent interactions of those particles are all subject to random processes. Thus, even if the primary particles are the same and have the same primary energies E₀ and the same zenith angles of incidence to the atmosphere, their shower development through the atmosphere will be different.

But it is the depth of the shower maximum, X_max, and the total number of electrons, N(X), at various depths X that will be different. Continuing our and other authors' works (see Section 2.1), we show that with respect to the shower maximum depth, the function describing fully the distributions of all electron characteristics is the same in any shower, independent of the primary energy, mass, angle of incidence, or fluctuations of shower development.

By referring to the shower maximum, we mean that instead of the depth X, one should use the age s of the shower at this depth, defined as

$$\begin{eqnarray}&&s(X)=3X/(X+2{X}_{\max }).\end{eqnarray} \tag{ 1 }$$

This was introduced by Hillas (1982) to describe a shower development stage analogously to a pure electromagnetic cascade where this particular relation comes out of the analytical solutions. It can be seen that the age s actually depends on the ratio $X/{X}_{\max }$, so that this very ratio could serve as a description of the shower development stage just as well. But we shall stick with tradition and use formula (1) for it.

The state of an electron at some age s can be described by four variables: the electron energy E, its distance to the shower axis in g cm⁻² or in $r/{r}_{M}$, where r_M is the Molière radius at the depth X(s), and the two angles of the electron direction with respect to the shower axis, the polar angle θ and the azimuth angle φ. Thus, the main conclusion we reach in this paper for the first time means that the five-dimensional function $f(\theta ,\varphi ,r/{r}_{M},E;s)$, representing the shape of the four-dimensional electron distribution at any level s, is the same for any shower, i.e., it has a universal dependence on the five variables (Section 2).

This conclusion could be drawn for showers having enough electrons in small cells of the four-dimensional space of the variables $E,r/{r}_{M},\varphi ,\theta$ at various levels s, so that the function $f(\theta ,\varphi ,r/{r}_{M},E;s)$ could be determined. This limits the shower primary energies to be greater than $({10}^{16}\mbox{--}{10}^{17})\,\mathrm{eV}$, depending on the choice of the variable cell size, and the range of the variable values to well-described ones. We show (Section 3) that some electron distributions can be described by a number of variables smaller than five, with the new variables being a product of the old ones raised to some power. This greatly simplifies the analytical parametrization of the electron distributions found here. In Section 4, we parametrize the simulated distributions by the analytical functions and give an estimation of the accuracy of our parametrization. In Section 5, we describe an analytical way of estimating the distortions of the electron angular and lateral distributions caused by the geomagnetic field.

The universality of the electron distributions in showers provides a new method for reconstructing shower development in the atmosphere using the fluorescence technique (Section 6). The full description of the electron states and its universal character allows one to predict exactly the light fluxes (both fluorescence and Cherenkov) emitted by shower electrons at consequent time intervals (as recorded by telescopes) by adopting only two quantities for a shower: N_max and X_max. Finding their values that best fit the data, one can determine quite well the shower primary energy, and from the depth X_max obtain information about the primary mass (see, e.g., Abraham et al. 2010a, 2010b; Abreu et al. 2013). A discussion and summary constitute Section 7.

2.1. Universality in Earlier Papers

2.2. General Considerations

The state of electrons in a shower is uniquely determined by the numbers ${\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E,s)$ given for all possibly occupied regions of the five variables. ${\rm{\Delta }}N$ is the number of electrons with angles to the shower axis $(\theta ,\theta +{\rm{\Delta }}\theta )$, with azimuth angle $(\varphi ,\varphi +{\rm{\Delta }}\varphi )$ around the axis parallel to that of the shower, drawn at the position of the electron which is at a lateral distance to the axis ($r/{r}_{M},r/{r}_{M}+{\rm{\Delta }}r/{r}_{M}$), with energy ($E,E+{\rm{\Delta }}E$), and at a shower age s. The Molière radius at a given height in the atmosphere is defined as ${r}_{M}\,=(21\,\mathrm{MeV}/\epsilon ){l}_{0}$, which is the mean square scattering angle of an electron with critical energy (∼82 MeV) along one radiation unit l₀. We define a function $f(\theta ,\varphi ,r/{r}_{M},E;s)$ as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E;s)=N(s)f(\theta ,\varphi ,r/{r}_{M},E;s)\\ &&\quad \times \,{\rm{\Delta }}\theta {\rm{\Delta }}\varphi {\rm{\Delta }}\mathrm{log}(r/{r}_{M}){\rm{\Delta }}\mathrm{log}E,\end{eqnarray} \tag{ 2 }$$

where N(s) is the total number of electrons at level s. Note that the function $f(\theta ,\varphi ,r/{r}_{M},E;s)$ has been defined for a single shower. However, we will show that for large showers, i.e., the number ${\rm{\Delta }}N$ is large ($\gt 100$), this function is universal. This means that it is independent of the primary particle energy, its mass, the incidence angle, or even the fluctuations in shower development (thus the requirement ${\rm{\Delta }}N\gt 100$). The function f can be represented as a product of functions, each depending on one variable and several parameters, as follows:

$$\begin{eqnarray}&&f(\theta ,\varphi ,r/{r}_{M},E;s)={f}_{E}(E;s){f}_{r}(r/{r}_{M};E,s){f}_{\varphi }\\ &&\quad \times \,(\varphi ;r/{r}_{M},E,s){f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s).\end{eqnarray} \tag{ 3 }$$

The variables on the right-hand side of the semicolon are the parameters of the function, whereas the independent variable is on the left side of the semicolon. The functions ${f}_{E},{f}_{r},{f}_{\varphi }$ and ${f}_{\theta }$ are normalized to unity when integrated over $\mathrm{log}\,E$, $\mathrm{log}(r/{r}_{M})$, φ, and θ respectively. They are, correspondingly, the distributions of the electron energy E at a given s, of the distance $r/{r}_{M}$ at given E and s, of the angle φ at given $r/{r}_{M}$, E, and s, and of θ at given φ, $r/{r}_{M}$, E, and s. As has been described above, the universality of the electron energy spectra and lateral distributions for various energies was already proven. It is seen from Equation (3) that to prove the universality of any electron distribution, represented by $f(\theta ,\varphi ,r/{r}_{M},E;s)$, it is left to demonstrate that each of the functions ${f}_{\varphi }(\varphi ;r/{r}_{M},E,s)$ and ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$ are universal.

To obtain the electron distributions, we simulated showers with CORSIKA (Heck et al. 1998), version 7.4, with the QGSJET-II model for high-energy interactions. For the atmospheric model, the US Standard Atmosphere was used. However, as we will later show, our conclusions do not depend on the above choices.

CORSIKA is a very elaborate computer code for simulating the development of air showers in the atmosphere. After choosing the primary particle, its energy, and incidence angle in the atmosphere, it follows all of the produced secondary particles (hadrons, muons, electrons, photons, and neutrinos) and their interactions. Of course, for hadrons, a high-energy interaction model has to be adopted. The most abundant particles in a shower are electrons (both signs), and it is their state, i.e., distributions over several variables at various atmospheric levels, that we are concerned with in this work.

We find the number of electrons ${\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E;s)$ in variable bins ${\rm{\Delta }}\theta {\rm{\Delta }}\varphi {\rm{\Delta }}(\mathrm{log}\,r/{r}_{M}){\rm{\Delta }}\mathrm{log}\,E$ at various ages s. Our study is restricted to variable regions that have the most electrons in the shower, i.e., to $0.7\lt s\lt 1.3$, 20 MeV < E < 200 MeV, and $r/{r}_{M}\lt 1.5$. The lower energy limit was chosen mainly based on our interest in describing the shower Cherenkov radiation with the threshold E = 21 MeV at sea level and increasing with height. For E < 20 MeV, accurate electron angular distributions are not necessary for the purpose of reconstructing the shower properties from the optical data because the fluorescence light is emitted isotropically. However, our parametrizations fit well down to ∼15 MeV.

Single showers with ${E}_{0}={10}^{17}\,\mathrm{eV}$ and 10¹⁶ eV were fully simulated, i.e., without using the thinning procedure (Hillas 1997), whereas for those with 10¹⁸ and 10¹⁹ eV, thinning (a procedure constraining the number of followed particles by assigning a weight to some of them) was, of course, used.

2.3. Distributions of Azimuth Angles

Figures 1 and 2 demonstrate the universality of the distributions ${f}_{\varphi }(\varphi ;r/{r}_{M},E,s)$ of the azimuth angle φ of electrons with two fixed energies E = 22 MeV and 220 MeV, at two distances $r/{r}_{M}$ for each E, at three shower ages $s=1,0.7$, and 1.3. The distributions refer to one proton shower with ${E}_{0}={10}^{19}\,\mathrm{eV}$ and to the average of 10 iron showers with ${E}_{0}={10}^{17}\,\mathrm{eV}$ (actually, the fluctuations from shower to shower for ${E}_{0}={10}^{17}\,\mathrm{eV}$ are not large, but for better comparison with the proton shower we have averaged them over 10 showers). We see that essentially there is no difference between the two curves in any of the graphs of Figures 1 and 2. The particular values of E and $\mathrm{log}(r/{r}_{M})$ have been chosen in such a way as to correspond to two values on both sides of the maximum of the $\mathrm{log}E$, and correspondingly $\mathrm{log}(r/{r}_{M})$, distributions (the latter describes the fraction of electrons in rings with a constant thickness ${\rm{\Delta }}\mathrm{log}(r/{r}_{M})=0.1$). At the maxima of these distributions, the independence of the primary particle characteristics is even better.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.4. Distributions of Polar Angles

Next, we move on to the next function, ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$, to be checked for universality—the distribution of the angle θ for fixed values of the parameters $\varphi ,r/{r}_{M},E,s$. Figures 3 and 4 illustrate the independence of the polar angle θ distributions of the electrons, ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$, from the primary particle energy and mass for s = 1. The chosen regions of the azimuth angle correspond to the electrons deflected mostly away from the shower axis, $0\lt \varphi \lt 20^\circ$; toward it, $160^\circ \lt \varphi \lt 180^\circ$; and perpendicularly to ${\boldsymbol{r}}$, $70^\circ \lt \varphi \lt 90^\circ$ and $\,90^\circ \lt \varphi \lt 110^\circ$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The average curve from 10 iron showers with the primary energy ${E}_{0}={10}^{16}\,\mathrm{eV}$ agrees very well with that referring to a single proton shower with 10¹⁹ eV. Comparing the average distribution with that for a single shower, one can also see fluctuations in individual bins, which are partly due to the small electron numbers in the bins (see Section 4.3). We have verified that such universality of the distributions ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$ holds at every stage of shower development in the range of $0.7\lt s\lt 1.3$. These angular distributions, when taken for the same values of energy and lateral distances of electrons (in units of the Molière radius), do show dependence on the shower development stage.

Studying the electron angular distributions in the multidimensional space of polar and azimuth angles, lateral distances, energies, and ages of showers, allows us also to search for a possible internal type of universality. This means that the angular distributions may depend on a combination of several of the five independent variables rather than on each of the variables separately. Such a universality would simplify the description of the angular distributions by reducing the number of necessary variables. A similar successful search describing the lateral distributions of electrons with various energies independent of their angles was done earlier (Giller et al. 2015).

3.1. Variable $\chi$

We found that introducing a new variable

$$\begin{eqnarray}&&\chi =(E/1\,\mathrm{MeV})\cdot r/{r}_{M}\end{eqnarray} \tag{ 4 }$$

the distributions ${F}_{\varphi }(\varphi ;\chi ,s)$, where ${F}_{\varphi }(\varphi ;\chi ,s)d\varphi$ gives the fraction of electrons with angles $(\varphi ,\varphi +d\varphi )$ at level s with $E\cdot r/{r}_{M}=\chi$, with respect to the total number of electrons with any φ, gives almost universal shapes. Thus, the function ${f}_{\varphi }(\varphi ;r/{r}_{M},E,s)$ reduces to a function ${F}_{\varphi }(\varphi ;\chi ,s)$ of three variables: $\chi ,\varphi$, and s. We have chosen E to be in MeV, so that χ is a dimensionless variable.

Figure 5 represents the universality of the distributions ${F}_{\varphi }(\varphi ;\chi ,s=1)$ of the azimuth angle φ of electrons with energies E = 22, 70, and 140 MeV and the corresponding lateral distances for two values of $\chi =7.8$ and 12.5. The curves are for one proton and one iron shower with energy 10¹⁷ eV without thinning. In Figure 6, we show how the universality holds in the case of a young (s = 0.7) and an old (s = 1.3) shower, this time taking electron distributions for one 10¹⁹ eV proton shower and ten 10¹⁷ eV iron showers with the thinning procedure switched on. The energies and lateral distances refer to $\chi =12.5$. Figures 5 and 6 also demonstrate the universal behavior of ${F}_{\varphi }(\varphi ;\chi ,s)$ with respect to the change of the primary particle energy or mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Variable $\xi$

It is not only the distributions ${F}_{\varphi }(\varphi ;\chi ,s)$ that are universal; the distributions of the polar angle θ also show a similar character. We have searched for one variable that could combine the energy, lateral distance, and θ angle for constant φ and s. We have found a variable ξ defined as

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{(E/1\mathrm{GeV})}^{\alpha }(1+r/{r}_{M})}{{(r/{r}_{M})}^{\beta }}\cdot \theta \ ,\end{eqnarray} \tag{ 5 }$$

where α and $\beta \gt 0$. We choose E to be in GeV and θ in degrees. Using this variable, we obtain for all regions of electron energies and lateral distances almost identical shapes of the normalized distributions of ξ for a given azimuth angle φ and shower age s. Thus, the function ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$ of five variables reduces to a function ${F}_{\xi }(\xi ;\varphi ,s)$ of only three variables: $\xi ,\varphi$, and s, where ${F}_{\xi }(\xi ;\varphi ,s)d\xi$ is the fraction of electrons with variable $(\xi ,\xi +d\xi )$ at angle φ at level s with respect to the total number of electrons with any ξ.

Using a procedure to minimize the differences between histograms binned for various energies and distances in different φ ranges, we have found the best values of α and β (see Equation (41) in the Appendix).

Figure 7 shows the distributions ${F}_{\xi }(\xi ;\varphi ,s)$ for two values of φ, at s = 0.7 and 1.3 (at s = 1, the picture is similar but with fewer fluctuations). For each φ, there are nine curves, referring to combinations of three electron energies, E = 23, 70, and 200 MeV, and three lateral distances, $\mathrm{log}(r/{r}_{M})=-1.25$, −0.85, and −0.45, for s = 0.7 and $\mathrm{log}(r/{r}_{M})=-1.05$, −0.65, and −0.25 for s = 1.3. The different sets of the lateral distances reflect the fact that in older showers, the lateral distribution becomes flatter. It is seen that the distributions are essentially independent of energy, lateral distance, or polar angle once the value of ξ is fixed.

Zoom In Zoom Out Reset image size

Having established the character of the universal distributions of the azimuth angles ${F}_{\varphi }(\varphi ;\chi ,s)$ and the distributions of the polar angles ${F}_{\xi }(\xi ;\varphi ,s)$, we undertake the task of finding their analytical descriptions.

4.1. Parametrization of the Distributions of the Azimuth Angles ${F}_{\varphi }(\varphi ;\chi ,s)$

The histograms in Figure 8 present examples of the distributions ${F}_{\varphi }(\varphi ;\chi ,s=1)$ in log-lin scale. The distributions look more or less like a fraction of a cosine function. Thus, we look for a function of the form

$$\begin{eqnarray}&&\mathrm{log}\,{F}_{\varphi }=A+B\cos (a\,\varphi +b),\end{eqnarray} \tag{ 6 }$$

where the parameters $A,B,a$, and b have to be fitted for different values of χ. We fit the distributions ${F}_{\varphi }(\varphi ;\chi ,s)$ for various χ in bins of ${\rm{\Delta }}\varphi =10^\circ$. The best-fitting polynomial descriptions, taking into account the dependence on s, are given in the Appendix.

Zoom In Zoom Out Reset image size

In Figure 8, the histogrammed distributions ${F}_{\varphi }(\varphi ;\chi ,s=1)$ refer to various values of χ for one 10¹⁷ eV iron shower. The smooth lines represent Equation (6) with the parameters described in Equations (38)–(39).

4.2. Parametrization of the Distributions of the Polar Angles ${F}_{\xi }(\xi ;\varphi ,s)$

We fit the distributions ${F}_{\xi }(\xi ;\varphi ,s)$ for various φ in bins of ${\rm{\Delta }}\varphi =10^\circ$ with the function

$$\begin{eqnarray}&&{F}_{\xi }(\xi ;\varphi ,s)=\displaystyle \frac{C\xi }{{(d+{e}^{\xi /2})}^{\gamma }}\,.\end{eqnarray} \tag{ 7 }$$

The parameters d and γ depend on φ, as is obvious from Figure 7. Parameter C normalizes the integral $\int {F}_{\xi }(\xi ;\varphi ,s)d\xi$ to unity. The best-fitting polynomial descriptions are in the Appendix.

Figure 9 presents an example of the distribution ${F}_{\xi }(\xi ;\varphi ,s)$ and fitting functions for four values of azimuth angles φ and three values of s. It can be seen that the description given by Equation (7) is quite satisfactory (see also Section 4.3).

Zoom In Zoom Out Reset image size

4.3. Estimation of the Accuracy of the Universal Fit

4.3.1. General Considerations

The analytical description of the electron distributions found in this paper, the universal fit, is not, of course, ideal. Now we aim to estimate how well it describes the actual electron distributions in an arbitrary shower detected by the fluorescence method as is done, e.g., at The Pierre Auger Observatory and other similar experiments.

The reasons why the actual distributions and the universal fit may differ are several:

• 1.  
  The universal fit is meant to describe an average from several simulated proton- and iron-initiated showers with primary parameters somewhat arbitrarily chosen, whereas a measured shower has a single primary energy E₀, primary mass, and zenith angle arising from the distributions of these parameters in the sample registered by the experiment.
• 2.  
  The simulated distributions have been fitted with possibly simple analytical functions, which might not be good enough in all regions of the variables.
• 3.  
  There is some inaccuracy in determining a particular shower age from the simulated shower, so that the derived electron distributions may refer to a slightly different age than the chosen one.
• 4.  
  Even for the same primary characteristics of the shower, the electron distributions may differ a bit above Poisson fluctuations in different showers due to some intrinsic fluctuations (see Section 4.3.3), although exact universality would mean that they do not.

4.3.2. Applicability to Auger Showers

First, we will check how well our universal fit describes the electron distributions in real (simulated) showers with primary energy ${E}_{0}={10}^{17}\,\mathrm{eV}$ at three ages s = 0.7, 1, and 1.3. To do this, we must build a sample to play the role of registered showers to be compared with our fit. We will no longer use the same sample used to derive the universal fit; it contains a smaller number of showers than we demand now, and they are all vertical and do not correspond to those registered by Auger. Treating it rigorously, we should simulate showers from the distributions of the primary energy E₀ and the zenith angle registered by Auger, with some adopted primary masses. However, a randomly chosen E₀ might be larger than the maximum energy in order to fully simulate a shower without the thinning procedure, which introduces rather large unwanted fluctuations. Moreover, our fit has been done for vertical showers, and although the electron component at the same atmospheric depth should not depend on the zenith angle, some small differences may be present due to, e.g., electrons from muon decay. Also, the primary masses are not known with good confidence.

Taking all of these into account, we constrain ourselves to showers with 10 primary energies drawn from the Gaussian distribution with $\sigma ={10}^{17}\times 0.14\,\mathrm{eV}$ (assuming the uncertainty in the energy scale to be 14%; Aab et al. 2015), the mean ${E}_{0}={10}^{17}\,\mathrm{eV}$, and some arbitrarily adopted primary mass composition: half protons and half iron nuclei. Showers with each energy are simulated at three zenith angles, 35°, 45°, and 55° (we checked that the Auger uncertainty of the zenith angle, ∼1°, affects distributions much less than E₀), so that the total number of showers in our sample equals n = 30.

For a fixed age s, we choose regions of the variables $(\varphi ,r/{r}_{M},E)$ where we want to check the accuracy of the fit of the θ distributions. These regions are chosen so as to contain relatively many electrons. We simulate the showers with the above primary characteristics and obtain the distributions ${\rm{\Delta }}{N}_{i}(\theta ;\varphi ,r/{r}_{M},E,s)$ for $i=1,2,\,\ldots .,\,n$.

To estimate how the distributions ${\rm{\Delta }}{N}_{i}(\theta ;\varphi ,r/{r}_{M},E,s)$ are described by the universal parametrized fit ${\rm{\Delta }}{N}_{\mathrm{par}}(\theta ;\varphi ,r/{r}_{M},E,s)$, we apply the ${\chi }^{2}$ test, the meaning of which in our case is the following.

The sample of 30 showers is believed to be drawn from the true distribution of ${\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E,s)$ (for each set of variables), with the true mean value

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{\mathrm{sh}}(\theta ,\varphi ,r/{r}_{M},E,s)\approx \langle {\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E,s)\rangle \end{eqnarray} \tag{ 8 }$$

being close to the mean from the sample, and the variance

$$\begin{eqnarray}&&{\sigma }_{\mathrm{sh}}^{2}(\theta ,\varphi ,r{/r}_{M},E,s)=\displaystyle \sum _{i=1}^{n}{({\rm{\Delta }}{N}_{i}-\langle {\rm{\Delta }}N\rangle )}^{2}/(n-1).\end{eqnarray} \tag{ 9 }$$

We want to check whether the parametrized values ${\rm{\Delta }}{N}_{\mathrm{par}}(\theta ,\varphi ,r/{r}_{M},E,s)$ can be described as drawn from the above distribution. In other words, here we treat ${\rm{\Delta }}{N}_{\mathrm{par}}(\theta ,\varphi ,r/{r}_{M},E,s)$ in the same way as a random number ${\rm{\Delta }}{N}_{i}(\theta ,\varphi ,r/{r}_{M},E,s)$.

Assuming that ${\rm{\Delta }}N(\theta ,\varphi ,r/{r}_{M},E,s)$ has a Gaussian distribution, the random variable ${\chi }^{2}$ defined as

$$\begin{eqnarray}&&{\chi }^{2}=\sum {({\rm{\Delta }}{N}_{\mathrm{par}}-{\rm{\Delta }}{N}_{\mathrm{sh}})}^{2}/{\sigma }_{\mathrm{sh}}^{2}(\theta ,\varphi ,r/{r}_{M},E,s),\end{eqnarray} \tag{ 10 }$$

where the sum is over all m cells of the variables $\theta ,\varphi ,r/{r}_{M},E,s$, has the ${\chi }^{2}$ distribution with $m-3$ degrees of freedom (two parameters, the mean and the variance, are determined from the sample data), provided that all m random values ${\rm{\Delta }}{N}_{\mathrm{par}}(\theta ,\varphi ,r/{r}_{M},E,s)$ are independent of each other (which, of course, is not quite true). We have taken into account only cells with ${\rm{\Delta }}N\gt 100$.

The ${\chi }_{1}^{2}$ for one degree of freedom equals

$$\begin{eqnarray}&&{{\chi }^{2}}_{1}=\displaystyle \frac{1}{m-3}\sum {({\rm{\Delta }}{N}_{\mathrm{par}}-{\rm{\Delta }}{N}_{\mathrm{sh}})}^{2}{/\sigma }_{\mathrm{sh}}^{2}.\end{eqnarray} \tag{ 11 }$$

We obtain ${{\chi }^{2}}_{1}=0.47$, meaning that the deflections of our fit are well-contained within the fluctuations between the chosen showers. So, the accuracy of the fit looks good enough for the reconstruction of the Auger showers at ∼10¹⁷ eV. However, we should keep in mind that the assumption of the independence of the fit values in the neighboring cells is not valid. The independence might be between groups of several cells, which is equivalent to a smaller number of degrees of freedom, resulting in turn to a larger ${\chi }^{2}$.

4.3.3. Compatibility with the True Distribution

As a measure of the accuracy of our parametrization independent of any experiment, we will use the mean relative difference between the parametrization predictions and the "true" values of ${\rm{\Delta }}{N}_{i}$, the latter of which is obtained from shower simulations for a fixed primary energy. We obtain

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{| {\rm{\Delta }}{N}_{\mathrm{par}}-{\rm{\Delta }}{N}_{\mathrm{sh}}| }{{\rm{\Delta }}{N}_{\mathrm{sh}}}\right\rangle =0.059,\end{eqnarray} \tag{ 12 }$$

where ${\rm{\Delta }}{N}_{\mathrm{sh}}$ are the previously used average numbers from the set of 30 showers (Section 4.3.2), and the average is taken over all cells.

To compare it with the minimum fluctuations of the multidimensional electron densities, i.e., those in a sample of simulated showers with the same primary parameters, we have simulated 10 vertical iron showers with ${E}_{0}={10}^{17}\,\mathrm{eV}$ and calculated the mean relative dispersion of ${\rm{\Delta }}{N}_{i}$ averaged over the variable space. Our result for all three ages is

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\sigma }_{\mathrm{sh}}}{{\rm{\Delta }}{N}_{\mathrm{sh}}}\right\rangle =0.064,\end{eqnarray} \tag{ 13 }$$

with values of 0.073, 0.051, and 0.071 for s = 0.7, 1, and 1.3, respectively.

Thus, the accuracy of our parametrization is roughly the same as the level of minimum shower fluctuations, ∼6%. We have used only those cells with ${\rm{\Delta }}{N}_{\mathrm{sh}}\gt 100$. Cells with a lower number of electrons do not have much physical meaning but add a lot to the dispersions ${\sigma }_{\mathrm{sh}}$ (Figure 10).

Zoom In Zoom Out Reset image size

The Poisson fluctuations in the chosen cells are less than 0.1. Subtracting them in each cell from ${\sigma }_{\mathrm{sh}}/{\rm{\Delta }}{N}_{\mathrm{sh}}$ in quadrature and averaging over all cells with ${\rm{\Delta }}{N}_{\mathrm{sh}}\gt 100$, we obtain for the intrinsic relative fluctuations

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\sigma }_{\mathrm{intr}}}{{\rm{\Delta }}{N}_{\mathrm{sh}}}\right\rangle =0.028,\end{eqnarray} \tag{ 14 }$$

with values of 0.038, 0.018, and 0.030 for s = 0.7, 1, and 1.3, respectively. Since the Poisson fluctuations have been subtracted, we expect that these numbers should not depend on the choice of the cell volume ${\rm{\Delta }}V$ as long as this volume is small. For proton-initiated showers, the fluctuations should be larger. For another choice of variable volume, one can add the corresponding Poisson fluctuation, obtaining the best accuracy with which fitting the real distributions makes sense.

The shower simulations we have done were with the geomagnetic field switched off, B = 0. With this field switched on, the directions of electrons would change and a new variable to describe the angular distributions of electrons would be added (Homola et al. 2015). As a result, the lateral distributions are also changed. The effect is not universal since it depends on the local magnetic field. It can be allowed for, if needed, by deforming the universal distributions. We will estimate how big this deformation is on particular levels in the atmosphere without shower simulations using a simple analytical method with some approximate assumptions. We will also give a recipe for obtaining corrected distributions from those without a ${\boldsymbol{B}}$ field when the effect is small.

Thus far, when we referred to electrons we have meant particles of both signs. In this section, however, "electrons" means negative-sign particles. For a deflection angle $\alpha \ll 1$, the new electron angles $\theta ^{\prime}$ and $\varphi ^{\prime}$ are

$$\begin{eqnarray}\theta ^{\prime} & = & \theta +{\rm{\Delta }}\theta =\theta +\cos \varphi \,\alpha ,\\ \varphi ^{\prime} & = & \varphi +{\rm{\Delta }}\varphi =\varphi -\displaystyle \frac{\sin \,\varphi }{\sin \theta }\,\alpha ,\end{eqnarray} \tag{ 15 }$$

with $\alpha \gt 0$ for ${e}^{-}$ and $\alpha \lt 0$ for ${e}^{+}$, assuming that $\cos \,\theta \sim 1$ and the azimuth angle φ is measured from the direction perpendicular to ${{\boldsymbol{B}}}_{\perp }$ clockwise. Thus, the new angular distribution functions ${f}_{\theta }^{{\prime} }$ and ${f}_{\varphi }^{{\prime} }$ are

$$\begin{eqnarray}{f}_{\theta }^{{\prime} }(\theta ;\varphi ,r/{r}_{M},E,s) & = & {f}_{\theta }[(\theta -\cos \,\varphi \,\alpha );\varphi ,r/{r}_{M},E,s]\,\mathrm{and}\\ {f}_{\varphi }^{{\prime} }(\varphi ;r/{r}_{M},E,s) & = & {f}_{\varphi }[\varphi +\displaystyle \frac{\sin \,\varphi }{\sin \theta }\,\alpha ;r{/r}_{M},E,s],\end{eqnarray} \tag{ 16 }$$

since ${\rm{\Delta }}\theta ^{\prime} ={\rm{\Delta }}\theta$ and ${\rm{\Delta }}\varphi ^{\prime} ={\rm{\Delta }}\varphi$.

A similar method is used to find the new lateral function ${f}_{r}^{{\prime} }$. If the new radial vector ${\boldsymbol{r}}^{\prime} ={\boldsymbol{r}}+{\boldsymbol{\rho }}$, then the change of the radial distance equals ${\rm{\Delta }}r=\rho \cos {\varphi }_{S}$, where $\rho \gt 0$ for ${e}^{-}$ and $\rho \lt 0$ for ${e}^{+}$, and ${\varphi }_{S}$ is the azimuth angle of ${\boldsymbol{r}}$ around the shower axis, measured from the direction perpendicular to ${{\boldsymbol{B}}}_{\perp }$ clockwise. Thus, the new function ${f}_{r}^{{\prime} }$ depends on the azimuth ${\varphi }_{S}$. Since all electrons with energy E at $r-\rho \cos {\varphi }_{S}$ on a surface with area ${\rm{\Delta }}S$ are shifted by the same distance ρ perpendicularly to ${{\boldsymbol{B}}}_{\perp }$, we have

$$\begin{eqnarray}&&\displaystyle \frac{{f}_{r}^{{\prime} }(r/{r}_{M},{\varphi }_{S};E,s)}{2\pi (r/{r}_{M})}{\rm{\Delta }}S^{\prime} \\ &&\quad =\,\displaystyle \frac{{f}_{r}(r/{r}_{M}-\rho /{r}_{M}\,\cos {\varphi }_{S};E,s)}{2\pi (r/{r}_{M}-\rho /{r}_{M}\,\cos {\varphi }_{S})}{\rm{\Delta }}S,\\ &&\quad \mathrm{where}\ {\rm{\Delta }}S^{\prime} ={\rm{\Delta }}S.\end{eqnarray} \tag{ 17 }$$

Thus, our recipe allowing for the geomagnetic effect in the lateral distribution function is the following:

$$\begin{eqnarray}&&{f}_{r}^{{\prime} }(r/{r}_{M},{\varphi }_{S};E,s)\simeq {f}_{r}(r/{r}_{M}\\ &&\quad -\,\rho {/r}_{M}\,\cos {\varphi }_{S};E,s)\left(1+\displaystyle \frac{\rho \,\cos {\varphi }_{S}}{r}\right).\end{eqnarray} \tag{ 18 }$$

Concerning the new functions for electrons of both signs, we have

$$\begin{eqnarray}{f}_{v}^{{\prime} }(v;\,...) & = & {P}^{-}(E;s){f}_{v}^{-^{\prime} }(v;\,...)\\ & & +{P}^{+}(E;s){f}_{v}^{+^{\prime} }(v;\,...),\end{eqnarray} \tag{ 19 }$$

where v stands for the variable r or $r/{r}_{M}$, φ, or θ; ${P}^{-/+}(E;s)$ are the fractions of electrons (positrons) with energy E at shower age s; and ${f}_{v}^{-^{\prime} }(v;\,...),{f}_{v}^{+^{\prime} }(v;\,...)$ are the above-determined functions for electrons $(\alpha ,\rho \gt 0)$ and positrons $(\alpha ,\rho \lt 0)$, respectively. At shower maximum, ${P}^{-/+}(E;s)$ are equal to 0.5 only at $E\gt 300\,\mathrm{MeV}$, with the fraction of positrons going down to $\sim 1/3$ at 10 MeV (see, e.g., Lafebre et al. 2009).

The values of the angular and lateral deflections α and ρ remains to be determined.

We will consider electrons and positrons with a given energy E on a level with age s. However, as has been already established (Giller et al. 2005a), the angular distributions of electrons of both signs with fixed energy do not depend on shower age s. They depend solely on the electron energy E, so that the shower age becomes irrelevant.

Hillas (1982) considered this problem by simulating low-energy electromagnetic cascades developing in the magnetic field, using detailed formulae for relevant cross-sections. For any electron energy E, he obtained an effective path length along which the ${\boldsymbol{B}}$ field should act, as if the electron had a constant energy E, to derive the correct mean deflection angle. For $E\gt 20\,\mathrm{MeV}$, the effective path lengths were a fraction of one radiation unit (r. u.), the largest being 0.73 r. u. at $E\to \infty$. Homola et al. (2015) studied the azimuthal asymmetry caused by the geomagnetic field of the electron angular distributions, integrated over the lateral spread and energy, in hadronic showers of the highest energies. They noted that the effect may be quite large in very inclined showers developing high in the atmosphere. Here, we study this effect for fixed lateral distances and fixed electron energies, where the angular and lateral distributions may be distorted, as it will turn out, by a rather small amount (which, however, may be necessary to take into account).

Considering an electron on its way back in a shower, one comes at some point to a photon, then to an electron (or positron), then again to an electron or positron, and so on. Photons are not deflected and positron deflection decreases the deflection accumulated by the electron. Calculating analytically the total deflection in this case, while taking into account the energy losses (even if possible), is not the aim of this paper. However, we can find the upper limit for the total deflection by considering the magnetic deflection of an electron maintaining its identity the entire time and losing energy via bremsstrahlung and ionization in an average way:

$$\begin{eqnarray}&&-d\,E/d\,t=E+\epsilon ,\end{eqnarray} \tag{ 20 }$$

with t being the path length in radiation units and $\epsilon =82\,\mathrm{MeV}$ the critical energy of the air. The solution is

$$\begin{eqnarray}&&E(l;{E}_{0})=({E}_{0}+\epsilon )\exp (-l/{l}_{0})-\epsilon \,,\,{\rm{}}\end{eqnarray} \tag{ 21 }$$

with $l/{l}_{0}=t$, where l₀ is the radiation unit in meters and E₀ the initial electron energy.

The magnetic deflection angle ${\rm{\Delta }}\varphi$ of a relativistic electron with energy E along the path length ${\rm{\Delta }}l\ll R(E)$, with R(E) being the curvature radius, equals

$$\begin{eqnarray}{\rm{\Delta }}\varphi & = & 0.03\displaystyle \frac{{B}_{\perp }({\rm{G}})}{E(\mathrm{MeV})}{\rm{\Delta }}l({\rm{m}})\\ & = & 1\buildrel{\circ}\over{.} 72\displaystyle \frac{{B}_{\perp }({\rm{G}})}{E(\mathrm{MeV})}{\rm{\Delta }}l({\rm{m}}).\end{eqnarray} \tag{ 22 }$$

Thus, on the way from E₀ to E, in a constant magnetic field, the electron is deflected by

$$\begin{eqnarray}\varphi & = & {B}_{\perp }\displaystyle \int \displaystyle \frac{{dl}}{E^{\prime} (l;{E}_{0})}\\ & = & {B}_{\perp }\displaystyle \frac{{l}_{0}}{\epsilon }\left(\mathrm{ln}\displaystyle \frac{E+\epsilon }{E}-\mathrm{ln}\displaystyle \frac{{E}_{0}+\epsilon }{{E}_{0}}\right).\end{eqnarray} \tag{ 23 }$$

For ${E}_{0}\to \infty$,

$$\begin{eqnarray}&&\varphi \to {B}_{\perp }\displaystyle \frac{{l}_{0}}{\epsilon }\mathrm{ln}\displaystyle \frac{E+\epsilon }{E}\\ &&\quad =\,\varphi (\epsilon ,{l}_{0})\mathrm{ln}\displaystyle \frac{E+\epsilon }{E}\\ &&\quad =\,5\buildrel{\circ}\over{.} 84{B}_{\perp }({\rm{G}})\mathrm{ln}\displaystyle \frac{E+\epsilon }{E}\,\end{eqnarray} \tag{ 24 }$$

at sea level, where $\varphi (\epsilon ,{l}_{0})$ denotes the deflection angle of an electron with the constant critical energy of the air along one radiation unit l₀, adopted as $36.2\,{\rm{g}}\,{\mathrm{cm}}^{-2}$.

Hillas determined the effective path length ${x}_{M}(E)$ in low-energy electromagnetic cascades, but due to the universality of electron distributions, the same path lengths should also apply to the high-energy hadronic showers. We will use them below.

The effective path length ${x}_{M}(E)$ as defined by Hillas equals, in our case,

$$\begin{eqnarray}{x}_{M}(E) & = & E\displaystyle \int \displaystyle \frac{{dl}}{E^{\prime} (l;{E}_{0})}\\ & = & {l}_{0}\displaystyle \frac{E}{\epsilon }\,\mathrm{ln}\displaystyle \frac{E+\epsilon }{E}\mathop{\to }\limits_{E\to \infty }{l}_{0}\,.\end{eqnarray} \tag{ 25 }$$

As expected, it is larger than that obtained by the cascade simulations of Hillas where ${x}_{M}(E)\mathop{\to }\limits_{E\to \infty }0.73\,{l}_{0}$, where the negatively charged electron progenitors are also photons and positrons not increasing their deflections.

The geomagnetic field at Earth is a fraction of 1 G. At the Auger site, B = 0.246 G. Due to its rather large declination angle of $-35^\circ$, the average B_⊥ for showers arriving isotropically with zenith angles smaller than 60° is not much smaller, ${B}_{\perp }\approx 0.9\,B\approx 0.22\,{\rm{G}}$. Thus, our estimation of the deflection upper limit for, e.g., E = 50 MeV would be $1\buildrel{\circ}\over{.} 6$ (adopting the air density at Auger of $1\cdot {10}^{-3}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$). The actual deflection, however, is smaller, and according to Hillas, the effective path length ${x}_{M}(E=50\,\mathrm{MeV})=12.35\,{\mathrm{gcm}}^{-2}$, which is a factor of 1.73 smaller than the upper limit (Equation (25)). Thus, the actual average deflection is $\alpha \sim 1^\circ$ at the Auger level. At higher atmospheric levels, $\alpha =1^\circ \exp ({h}_{A}/8.5\,\mathrm{km})$, where h_A is the height above the Auger site ($\sim 1.4\,\mathrm{km}\,{\rm{a}}.{\rm{s}}.{\rm{l}}.$).

At the same time, the electron would acquire via Coulomb scattering an angle θ with respect to the shower axis, which we can estimate by a similar way to that above. Adopting the average increase rate of $\langle {\theta }^{2}\rangle$ as

$$\begin{eqnarray}&&\displaystyle \frac{\langle d\,{\theta }^{2}\rangle }{{dt}}={\left(\displaystyle \frac{21}{E(\mathrm{MeV})}\right)}^{2}\,\,{\rm{}}\end{eqnarray} \tag{ 26 }$$

and the energy losses from Equation (20), we obtain for ${E}_{0}\to \infty$

$$\begin{eqnarray}&&\sqrt{\langle {\theta }^{2}\rangle }=\displaystyle \frac{21\,\mathrm{MeV}}{\epsilon }\sqrt{\displaystyle \frac{\epsilon }{E}-\mathrm{ln}\left(\displaystyle \frac{E+\epsilon }{E}\right)}\,,\end{eqnarray} \tag{ 27 }$$

which gives $\sqrt{\langle {\theta }^{2}\rangle }\sim 12^\circ$ for E = 50 MeV. Thus, in comparison with the above calculated $1\buildrel{\circ}\over{.} 6$ magnetic shift of such an electron, its Coulomb scattering looks quite big. It is worth noting, however, that for actual electrons in a shower, both angles will be smaller, but the Coulomb to magnetic ratio will increase since the former decreases due to the parent photons and the latter decreases due to photons and positrons.

From above, we can draw the conclusion that the effect is rather small for typical electron parameters, at least at heights ${h}_{A}\lt 20\,\mathrm{km}$ where $\alpha \lesssim 10^\circ$, so that we can find the relative change of the angular distribution as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}({f}_{\varphi }{f}_{\theta })}{{f}_{\varphi }{f}_{\theta }}=\displaystyle \frac{{\rm{\Delta }}{f}_{\varphi }}{{f}_{\varphi }}+\displaystyle \frac{{\rm{\Delta }}{f}_{\theta }}{{f}_{\theta }}\\ &&\quad =\,-\displaystyle \frac{1}{{f}_{\varphi }}\displaystyle \frac{\partial {f}_{\varphi }}{\partial \varphi }\displaystyle \frac{\sin \,\varphi }{\sin \,\theta }\alpha +\displaystyle \frac{1}{{f}_{\theta }}\displaystyle \frac{\partial {f}_{\theta }}{\partial \theta }\cos \,\varphi \,\alpha \,,\end{eqnarray} \tag{ 28 }$$

where we have assumed that the number of electrons ${\rm{\Delta }}N(r/{r}_{M};E,s)$ stays unchanged. The relative change of the angular distribution function depends strongly on the choice of angles where the derivatives are taken (see Figures 1–4). For some representative values, $s=1,E=50\,\mathrm{MeV},\mathrm{log}(r/{r}_{M})\,\sim -0.7$, and regions of angles where the derivatives are more or less constant, we get

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}({f}_{\varphi }{f}_{\theta })}{{f}_{\varphi }{f}_{\theta }}=\left(0.025\displaystyle \frac{\sin \,\varphi }{\sin \,\theta }-0.085\,\cos \,\varphi \right)\alpha ,\,{\rm{}}\end{eqnarray} \tag{ 29 }$$

where α is in degrees.

Using the effective path length for an electron with energy E, the lateral shift ρ equals

$$\begin{eqnarray}&&\rho =\displaystyle \frac{1}{2}\displaystyle \frac{{l}_{\mathrm{eff}}^{2}(E)}{R(E)}\,,\end{eqnarray} \tag{ 30 }$$

with R(E) being the electron curvature radius. Thus, an electron with E = 50 MeV and ${l}_{\mathrm{eff}}=12.35\,{\rm{g}}\ {\mathrm{cm}}^{-2}$ (Hillas 1982) will be deflected on average by $\rho \sim 1$ m at the Auger level. Using the shape of the lateral distribution function f_r (Giller et al. 2015), we obtain for s = 1 and $\mathrm{log}(r/{r}_{M})=-0.7$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{f}_{r}}{{f}_{r}}=0.003\,\cos \,{\varphi }_{S}\,\rho ({\rm{m}}){\rm{.}}\,{\rm{}}\end{eqnarray} \tag{ 31 }$$

The total relative change of the number of electrons for our example, adopting $\theta =10^\circ$, equals

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}({f}_{r}{f}_{\varphi }{f}_{\theta })}{{f}_{r}{f}_{\varphi }{f}_{\theta }}=0.003\cos {\varphi }_{S}\,\rho ({\rm{m}})\\ &&\quad +\,(0.144\sin \varphi -0.085\cos \varphi )\alpha (^\circ ),\end{eqnarray} \tag{ 32 }$$

where for the Auger site $\rho =\pm 1\,{\rm{m}}\,\exp ({h}_{A}/8.5\,\mathrm{km})$, $\alpha \,=\pm 1^\circ \,\exp ({h}_{A}/8.5\,\mathrm{km})$, and the + sign is for ${e}^{-}$. For electrons of both signs, we have

$$\begin{eqnarray}\displaystyle \frac{{\rm{\Delta }}({f}_{r}{f}_{\varphi }{f}_{\theta })}{{f}_{r}{f}_{\varphi }{f}_{\theta }} & = & {P}^{-}(E;s)[0.003\cos {\varphi }_{S}\,\rho ({\rm{m}})\\ & & +(0.144\sin \varphi -0.085\cos \varphi )\alpha (^\circ )]\\ & & +{P}^{+}(E;s)[-0.003\cos {\varphi }_{S}\,\rho ({\rm{m}})\\ & & -(0.144\sin \varphi -0.085\cos \varphi )\alpha (^\circ )]\,,\end{eqnarray} \tag{ 33 }$$

The two angles ${\varphi }_{S}$ and φ are, however, correlated. The average electron direction vector at a given distance ${\boldsymbol{r}}$ lies in the plane containing ${\boldsymbol{r}}$ and the shower axis, so that $\langle \varphi \rangle ={\varphi }_{s}$. Adopting both angles to be equal, we obtain for the maximum value of the total relative change

$$\begin{eqnarray}&&{\left|\displaystyle \frac{{\rm{\Delta }}f}{f}\right|}_{\max }\approx 0.033\exp ({h}_{A}/8.5\,\mathrm{km})\displaystyle \frac{{B}_{\perp }}{0.22\,{\rm{G}}}\,,\end{eqnarray} \tag{ 34 }$$

where $f=f(\theta =10^\circ$, $\varphi \lt 40^\circ$, $r/{r}_{M}={10}^{-0.7}$, $E=50\,\mathrm{MeV},s=1)$, ${P}^{-}\,(50\,\mathrm{MeV};s\,=\,1)=0.6$, and ${P}^{+}(50\,\mathrm{MeV};s\,=\,1)\,=0.4$. A more representative quantity is, however, the average modulus of the right-hand side of Equation (33) when the numerical coefficient equals $0.033\cdot 2/\pi =0.021$. Thus, for our example, the relative change of the total distribution function is below ∼20% up to heights ${h}_{A}\lesssim 20\,\mathrm{km}$.

The geomagnetic correction depends on the value of

$$\begin{eqnarray}&&\exp ({h}_{A}/8.5\,\mathrm{km})\displaystyle \frac{{B}_{\perp }}{0.22\,{\rm{G}}}\propto \displaystyle \frac{{B}_{\perp }}{{\rho }_{\mathrm{air}}(h)}=a,\end{eqnarray} \tag{ 35 }$$

where ${\rho }_{\mathrm{air}}(h)$ is the air density at height h of the shower level. The parameter a derived here is just the same as that used by Homola et al. (2015) for a parametrization of their simulation results. They calculate the asymmetry of the angular distribution of Cherenkov photons (being almost the same as that of electrons) integrated over lateral distance and electron energy. For example, for a 10¹⁹ eV shower maximum at heights 8–10 km with ${B}_{\perp }=0.2\,{\rm{G}}$, they obtain ${\left|\tfrac{{\rm{\Delta }}f}{f}\right|}_{\max }\approx 0.2$, whereas for our example parameters (Equation (34); allowing for a slightly different θ) we get ∼0.1. Our result, however, can be treated as a lower limit since we have chosen the variable regions where the distribution functions have the second derivative almost zero, so that if the numbers of electrons and positrons were equal the effect would be zero. Since the result of Homola et al. refers to an average effect for all electrons emitting Cherenkov light ($E\gt 30\,\mathrm{MeV}$ at that height), both results are compatible.

The universality of the electron distributions in showers opens a way to a new method for reconstructing shower development in the atmosphere, as achieved with the fluorescence technique (Abraham et al. 2010c). At the Pierre Auger Observatory, the fluorescence detector measurements consist of shower light images registered within small time bins by a (at least one) camera with 440 PMT pixels with a $1\buildrel{\circ}\over{.} 5$ field of view each. Shower reconstruction means determining the shower profile, i.e., the number of charged particles, N(X), as a function of the atmospheric depth X from the measured time series of the light images. Since the shape of N(X), having the form proposed by Gaisser and Hillas (Gaisser & Hillas 1977), changes little from shower to shower, it is N_max and X_max that determine it to a good approximation. Adopting their values and using the universal description of the electron states at any shower age allow one to predict exactly the light fluxes (both fluorescence and Cherenkov) emitted by shower electrons at consequent time intervals. We would like to stress that to do it, the full function $f(\theta ,\varphi ,r/{r}_{M},E;s)$ is needed for showers with measurable lateral size. Finding values of N_max and X_max that best fit the data, one can determine the shower primary energy and from the depth X_max obtain an estimate of the primary mass.

This method would be particularly useful for reconstructing showers with a small detecting angle (angle between shower axis and camera symmetry axis), since they contain a non-negligible amount of Cherenkov light as these are registered by an additional set of telescopes (at Auger), HEAT (Aab et al. 2015). It is thanks to finding the universal function $f(\theta ,\varphi ,r/{r}_{M},E;s)$ in this paper that this light can be treated accurately for the first time.

This paper is the last in a series where gradually more variables describing the electron state in a shower are taken into account. First, it was the electron energy spectra that were shown to be universal and dependent on shower age only (Giller et al. 2004). Then, the lateral distributions expressed in Molière radius were shown to be universal as well (Giller et al. 2005b; Nerling et al. 2006). Giller et al. (2005b) studied the lateral distributions for electrons with fixed energies and found that they are independent even of the shower age and could be described by a universal, one-dimensional function.

Finally, the angular distributions, dependent on the electron energy and lateral distance, have been found in this paper for the first time. As one would expect, they are universal as well. The universality has made the analytical description of electron distributions meaningful and relatively easy. In particular, the number of variables describing them has been found to be smaller than five, disclosing a new sort of universality as an inner characteristic of showers (Giller et al. 2015 called it an extended universality).

Universality here means independence from the primary particle energy or mass, zenith angle, or fluctuations from shower to shower. We do not mention the high-energy interaction model since the independence from the primary mass implies independence from the model. Indeed, one could imagine that the iron shower is actually a proton one, with the interaction model corresponding to the iron nucleus. This would be a dramatic change for the proton interactions, much bigger than the differences between the existing models. A demonstration of this could be the dependence of X_max at some primary energy on the interaction model and primary mass. The difference in X_max between proton and iron showers is considerably larger than that for different models and the same primary mass.

The universality of the angular distributions, dependent on shower age, electron energy, and lateral distance, has been checked in the primary energy region ${10}^{16}\mbox{--}{10}^{19}\,\mathrm{eV}$. However, since the lateral distributions are universal for energies up to 10²⁰ eV, it would be hard to believe that the angular distributions are not. For energies lower than 10¹⁶ eV, the numbers of electrons in the cells ${\rm{\Delta }}\mathrm{log}E\,{\rm{\Delta }}\mathrm{log}(r/{r}_{M}){\rm{\Delta }}\varphi {\rm{\Delta }}\theta$ would mostly be too small to determine the density there. Nevertheless, there seems to be no reason why universality should stop at 10¹⁶ eV, although it could be checked for averages of many showers only and not for single showers, as could be done at higher energies.

The distributions found here are normalized to unity. Thus, to find the actual number of electrons in a given cell of variables at a particular depth X in the atmosphere, one has to know the total number of electrons at this depth, N(X). Without going too much into details, this can be done by adopting the two free parameters of the Gaisser–Hillas function, N_max and X_max, describing N(X).

The universality of showers, consisting mainly of electrons, opens a new method for deriving N(X) from the optical data of observatories such as The Pierre Auger Observatory and Telescope Array. They measure the light fluxes while it is propagating through the atmosphere, and it is not only the fluorescence but also the Cherenkov light. From the comparison of the measured light intensities with those predicted, the best fitting N_max and X_max can be found. However, to predict well both light components, the full electron distribution $f(\theta ,\varphi ,r/{r}_{M},E,s)$ has to be known. It is particularly important for the HEAT addition to Auger, which measures closer, lower energy showers where the Cherenkov light plays a much bigger role than in the much distant Auger showers. The distributions found here, describing the electron angular distributions as function of their lateral distance and energy, allow one to predict exactly the Cherenkov contribution for the first time.

In some cases, a distortion of the angular and lateral distributions by the geomagnetic field would have to be taken into account. We have derived analytically when and how to allow for it.

Our parametrized distributions diverge from the simulated ones typically by ∼6%. It is roughly the same as the dispersion between showers with the same primary conditions. It is also well within the experimental uncertainties resulting from the accuracy of the determination of the primary energy at the Pierre Auger Observatory.

We are particularly grateful to the authors of the CORSIKA computer code, since without it our results would not have been achieved. We thank the Pierre Auger Collaboration who showed us the importance of this work. We also thank the referee whose suggestions have made the paper much improved. This work has been supported by grant No. DEC-2013/10/M/ST9/00062 of the Polish National Science Center.

A.1. Distributions of the Azimuth Angles ${F}_{\varphi }(\varphi ;\chi ,s)$

Using the variable

$$\begin{eqnarray}&&\chi =E\cdot r/{r}_{M},\end{eqnarray} \tag{ 36 }$$

where E is in MeV, we have reduced the distribution ${f}_{\varphi }(\varphi ;r/{r}_{M},E,s)$ of four variables to the distribution ${F}_{\varphi }(\varphi ;\chi ,s)$ of three variables, so that

$$\begin{eqnarray}&&{f}_{\varphi }(\varphi ;r/{r}_{M},E,s)={F}_{\varphi }(\varphi ;\chi =E\cdot r/{r}_{M},s).\end{eqnarray} \tag{ 37 }$$

We fit the distributions ${F}_{\varphi }$ with the function

$$\begin{eqnarray}&&\mathrm{log}\,{F}_{\varphi }=A+B\cos (a\,\varphi +b),\end{eqnarray} \tag{ 38 }$$

where φ is in degrees and the best fitted parameters are

$$\begin{eqnarray}&&\begin{array}{l}A\,=\,-2,\\ B\,=\,(0.491-0.084\cdot s)+(0.482-0.178\cdot s)\cdot \mathrm{log}\,\chi ,\\ a\,=\,(0.817-0.290\cdot s)+(\,-0.035+0.144\cdot s)\cdot \mathrm{log}\,\chi ,\\ b\,=\,(58.44+22.92\cdot s)+(0.814-11.80\cdot s)\cdot \mathrm{log}\,\chi .\end{array}\end{eqnarray} \tag{ 39 }$$

A.2. Distributions of the Polar Angles ${F}_{\xi }(\xi ;\varphi ,s)$

Using the variable

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{E}^{\alpha }(1+r/{r}_{M})}{{(r/{r}_{M})}^{\beta }}\cdot \theta ,\end{eqnarray} \tag{ 40 }$$

where E is in GeV, θ is in degrees, and the parameters α and β are

$$\begin{eqnarray}\alpha & = & 0.61+0.76\cdot {10}^{-3}\cdot \varphi \,{\rm{a}}{\rm{n}}{\rm{d}}\\ \beta & = & 0.46-2.1\cdot {10}^{-3}\cdot \varphi \\ & & +5.5\cdot {10}^{-6}\cdot {\varphi }^{2}\,{\rm{w}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{e}}\,\varphi \,{\rm{i}}{\rm{s}}\,{\rm{i}}{\rm{n}}\,{\rm{d}}{\rm{e}}{\rm{g}}{\rm{r}}{\rm{e}}{\rm{e}}{\rm{s}},\end{eqnarray} \tag{ 41 }$$

we have reduced the distribution ${f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)$ of five variables to the distribution ${F}_{\xi }(\xi ;\varphi ,s)$ of only three variables, so that

$$\begin{eqnarray}&&{f}_{\theta }(\theta ;\varphi ,r/{r}_{M},E,s)\\ &&\quad =\,{F}_{\xi }\left(\xi =\displaystyle \frac{{E}^{\alpha }(1+r/{r}_{M})\theta }{{(r/{r}_{M})}^{\beta }};\varphi ,s\right)\cdot \displaystyle \frac{d\,\xi }{d\,\theta }.\end{eqnarray} \tag{ 42 }$$

We fit the distributions ${F}_{\xi }(\xi ;\varphi ,s)$ with the function

$$\begin{eqnarray}&&{F}_{\xi }(\xi ;\varphi ,s)=\displaystyle \frac{C\xi }{{(d+{e}^{\xi /2})}^{\gamma }}\,,\end{eqnarray} \tag{ 43 }$$

where C is the normalization constant and the best fitted parameters are

$$\begin{eqnarray}\mathrm{log}\,d & = & {d}_{1}-{d}_{2}\cdot {10}^{-3}\cdot \varphi +{d}_{3}\cdot {10}^{-5}\cdot {\varphi }^{2}\ ,\\ \mathrm{where}\,{d}_{1} & = & 1.177-0.066\cdot s-0.178\cdot {s}^{2}\ ,\\ {d}_{2} & = & 15.3-0.82\cdot s-2.86\cdot {s}^{2}\ ,\\ {d}_{3} & = & 2.28-1.48\cdot s+0.31\cdot {s}^{2}\\ \mathrm{and}\,\gamma & = & 1.60+1.85\cdot {10}^{-2}\cdot \varphi -3.92\cdot {10}^{-5}\cdot {\varphi }^{2}\ .\end{eqnarray} \tag{ 44 }$$