1 Introduction

The basic method of extensive air shower registration is the use of surface arrays, usually consisting of scintillation detectors (KASCADE, HiRes, TA, AGASA, etc.), resistive plate chambers (e.g. ARGO-YBJ) and sometimes Cherenkov detectors (Pierre Auger Observatory, LHAASO, HAWC). More recently, small extensive air shower instruments have become popular, recording small showers initiated by particles with energies below the knee (\(E<10^{15}\) eV) (HiSPARC, QuarkNet, WALTA, NALTA, ALTA, SALTA, CZALTA, SKALTA, CHICOS, CROP, CosMO, \(\mu \)Net, EEE, ADA, GELATICA, \(\mu \)Cosmics, Astroneu, EAS-UAP, CREDO-Maze). Their applications are mostly for educational and popularisation purposes. Sometimes, however, they are built with the idea of searching for phenomena that have not yet been observed, but are expected or at least possible in theory. These may be correlated distant showers initiated by nuclei in the Gerasimova–Zatsepin process [1], or cosmic rays from extremely energetic massive hypothetical particles, or ultrahigh-energy photons of unknown origin.

Due to the very steep energy spectrum of the primary cosmic rays, an array of several small detectors, closely spaced, with a detection threshold of 1 m.i.p. for each detector, can mostly respond to very local shower reaching the ground with a small total number of particles, even in the hundreds and less. In the showers recorded by such instruments, the local expected particle densities at the detector locations are generally less than the set threshold. Triggers appear as a result of the upward fluctuation in the number of particles hitting the detector. Therefore, knowledge of the distribution of the number of particles hitting a detector with a given area at a given location is necessary to determine at least the count rate of the entire apparatus. At the same time, the determination of the density distributions in the detector alone is not sufficient to determine the counting rate. For the formation of the trigger of a small shower array, it is necessary to examine the correlations between the particles in the shower. For a complete picture of the mechanism of formation of a specific trigger, we should know the magnitude of the two-particle correlation as a function of the distance between the particles.

We are used to thinking of extensive air showers as events made up of large numbers of individual particles. Most shower arrays are cascades of millions or more particles, and we expect that, although they are genetically related, the correlations between them will be lost by the law of large numbers.

On the other hand, almost every high-energy electron recorded by the detector at the Earth’s surface comes as a piece from an electron/positron pair created by a high-energy photon in the process of pair creation. This process must have taken place not high above the detector, before the electron and positron interacted with atmospheric nuclei to produce further photons, and these further electrons would have blurred the original correlation, which for this reason alone should be short range (compared to the radial size of the entire shower).

Both effects, the apparent short-range correlation and its blurring by many other, already significantly less correlated particles, define what we might call the density fluctuations of the showering particles. For independent particles, the density fluctuations are of course Poissonian, and by adding up the total number of particles in the shower, they give the Poissonian fluctuations of Bhabha and Heitler [2]. In general, the distribution of the number of particles in a particular detector need not be so. Any correlation will result in deviations from the Poisson distribution, and these deviations are those that we will be looking for.

The dominant fluctuations we need to consider are those associated with the stochastic nature of the extensive air shower developing in the atmosphere. These are, first, fluctuations in the shower size, by which we refer to the total number of particles at the observation level, and fluctuations of the shape of the transverse distribution function of the particles in the shower. We assume the axial symmetry of the distribution function, and its form proposed by Greisen known as the Nishimura–Kamata–Greisen (NKG) function [3, 4].

These correlations say nothing about the mechanisms involved in the physical, generic relationship between the particles in the cascade, and we will not be interested in it. The appropriate method of analysis will circumvent this problem.

2 Method of the correlation analysis

The study of correlations is one of the key problems in high-energy physics in trying to describe multiparticle production processes in strong interactions. Hadronisation and multiplicative production processes in general resemble to a large extent the cascade process of shower development. Self-similar multiplicative branching models have led to many original methods for the analysis and description of self-similar fluctuations over a wide range of scale parameter variations.

The density of particles in a shower is described by the function \(\rho (\textbf{r})\), where the vector \(\textbf{r}\) determines the location and possibly other parameters of interest of the particle. Once normalised to the total size of the shower and restricted to the location of the particles in the observation plane, we can define the probability density of finding a particle at a given location q(x), where x is actually a two-dimensional vector.

Correlations arise when we look at the positions of several shower particles simultaneously.

If we are interested in the probability of finding a pair of particles in a detector of size \(\Delta \), we define a two-dimensional density \(\rho _2(x_1, x_2)\), and the quantities we will use to determine the interparticle correlations are normalised cell-averaged factorial moments. The scaled factorial moments were first introduced by Bialas and Peschanski [5, 6].

$$\begin{aligned} F_2(\Delta ) = \frac{\int _{\Omega s}\rho _2(x_1, x_2)dx_1 dx_2 }{ \int _{\Omega s} \rho (x_1)\rho (x_2) dx_1 dx_2 } \end{aligned}$$
(1)

where \(\Omega s\) is the subspace of the set of all possible 2-particle positions satisfying the condition that each of its points corresponds to the case where both particles hit a detector of fixed size \(\Delta \). In general, we can define the k-particle density \(q_k(x_1, x_2,\ldots , x_k)\) and similarly the corresponding factorial moments of higher orders.

The definition of the factorial moments given in Eq. (1) is very convenient for practical application. Determining the values of the corresponding integrals boils down to counting the pairs (threes, fours,...) of shower particles that satisfy the condition that the usual spatial distance between them is less than \(\Delta \) [7, 8].

$$\begin{aligned} F_k(\Delta )= & {} \frac{1}{\text {norm.}}\ \ k! \ \ \sum \limits _{i_1<i_2<\ldots <i_k} \ \ \nonumber \\{} & {} \quad \prod \limits _{ \begin{array}{c} {\scriptstyle \text {all pairs}} \\ \scriptstyle m_1, m_2 \end{array} }\normalsize \Theta \left( \Delta {-} {\text {distance}}\left( x_{i_{m_1}},x_{i_{m_2}}\right) \right) \end{aligned}$$
(2)

where \(\Theta \) is the Heaviside unit-step function.

Counting the particles in the shower obtained from the simulation program is, in general, not a problem. The main difficulty is the factor called the “norm” in equation (2). As defined in eq. (1), this is the same number of occurrences of pairs (threes, fours, etc.) that fit into the detector if all the shower particles that the detectors could potentially detect are independent of each other.

The construction of a reference sample of showers free from correlation is impossible by the very nature of the problem, due to the presence of non-Poissonian fluctuations in the size and distribution of the particles in the showers. As these two sources of particle correlation are dominant and are not the subject of this work, we have eliminated them, leaving both the shower size and the particle transverse distribution in the reference sample in each event the same as in the simulated case. The independence of the shower particles was ensured by the random rotation of each particle around the intersection of the shower axis and the observation plan according to the assumed azimuthal symmetry of the shower.

The reference shower formed in this way retains the correlations associated with age and radial distribution and, as the denominator in the formula Eq. (1), ensures that the determined values of the factorial moments \(F_k\) show only correlations due to other factors included in the evolution of the shower described by the shower simulation program.

3 Results

The simulated showers used to determine the expected effect of particle correlation in the shower were generated using the CORSIKA program [9]. It was originally developed over 30 years ago in Karlsruhe for the KASCADE experiment [10,11,12] and is still one of the most popular programs for simulating extensive air showers, even at the highest observed energies.

Among the many options available to the user in shower evolution calculations, one of the most important appears to be the choice of interaction models. In the high-energy region \(\sqrt{s} \ge \)100 GeV, we have basically used the EPOS LHC [13] (v3400) model. We also checked the results obtained with the SIBYLL 2.3d, VENUS and QGSJET-II-04 models and found no significant differences. This is not surprising, since our interest is focused on primary particle energies of 10\(^{11}\)–10\(^{13}\) eV, which were available in accelerator experiments more than a half century ago, and all high-energy interaction models implemented in CORSIKA were from the beginning matched to accelerator results at these energies. Rather, we expect that the interaction model used to simulate interactions at lower energies may have an impact on the results of the analysis. In this respect, CORSIKA is equipped with three options: GHEISHA, FLUKA (FLUktuierende KAskade) and UrQMD (Ultrarelativistic Quantum Molecular Dynamics) model.

It is reasonable that known electromagnetic processes are responsible for the local correlations of particles in bursts. These are described in CORSIKA by the EGS4 (Electron Gamma Shower) package [14].

Showers initiated by protons with energies up to 100 GeV mostly do not contain any charged particles at sea level. They extinguished before reaching the observation plane. Since we are also interested in cases where several charged particles are in close proximity, and the Molière unit describing the scale of longitudinal particle propagation is on the order of 100 ms, we had to run millions of shower simulations.

Based on these statistics, we determined the normalised factorial moments of order 2, 3 and 4 for reciprocal distances \(\Delta \) from a few tens of centimetres to a few tens of metres according to Eq. (2). Below this range, the statistics at the lowest energies were too small to obtain statistically significant results anyway, and the resulting significance for the results of possible and potential measurements is in any case marginal. For distances above a few tens of metres, the correlations between particles are no longer local, and the effects associated with the adopted method of creating a “mixed events” database begin to become apparent. Particles randomly rotated around the axis of the shower, maintaining the shape of the transverse distribution, remain in the majority at distances of tens of metres.

The results of these calculations are shown in Fig. 1.

Fig. 1
figure 1

Second (a), third (b) and fourth (c) normalised factorial cumulants for electron component in CORSIKA-simulated showers calculated for different initial particle proton energy. Thin straight lines along the thick simulation results are the power-law fits in the region 1 to 10 ms

Full size image

The methodology used is to measure the dependence of the normalised factor moments \(F_k\) on the size of the k-particle cluster \(\Delta \). A particle scattering that does not exhibit fluctuations other than statistical (Poisson) has the property that \(F_k(\Delta \)) is independent of the resolution of \(\Delta \) (in the limit of \(\Delta \rightarrow 0\)). On the other hand, if there is a correlation and it is self-similar (“intermittent”), which would be expected in the cascade regime, then \(F_k\) obeys a power law of

$$\begin{aligned} F_k(\Delta )= (\Delta )^{-\phi (k)}~~,~~~~~~(\Delta \rightarrow 0) \end{aligned}$$
(3)

However, the attempt to interpret extensive air showers as fractal objects and to analyse their topological dimensions on this basis is not entirely legitimate.

As shown in Fig. 1a, normalised factorial moments of order 2 show for all analysed energies a clear power-law character. In addition to the results obtained for the simulated CORSIKA showers shown by thicker lines in the figures, the power-law fits to these results are shown by (sometimes hard-to-see) straight thin lines.

The values of the power-law index \(\phi \) for the second and third moments are very close to the value of 1. For the second moment (Fig. 1a), they oscillate around 0.9. For the third moment (Fig. 1b), its values remain around 0.7. For the fourth moment (Fig. 1c), the values oscillate around 0.3.

4 Decoherence curve

Correlations among the secondary products of extensive air showers have been studied since their discovery by Auger [15, 16]. Since then, the coincidence counting rate of two cosmic ray detectors as a function of the separation between them has been measured under various conditions and up to large relative distances.

Fig. 2
figure 2

Normalised number of coincidences versus distance between detectors measured in [17] (black dots), [18] (empty circles) and [19] (squares). The line represents \(\Delta ^{-1}\)

Full size image

Figure 2 presents the decoherence curve recently measured by Riggi and co-workers [17] and, for comparison, the old historical measurements of Auger [18] and Kolhörster [19]. Riggi used an array of 15 telescopes consisting of two \(20\times 20\) cm\(^2\) scintillation counters placed one above the other and recorded the number of coincidences of each pair of telescopes for almost a month. The results show a clear power-law pattern of decreasing coincidence counts for distances up to about 5 m distance between the detectors. The slope of the dependency plotted in Fig. 2 by the solid line is equal to 1. The different experimental results were obtained under different conditions, at different elevations above sea level and with different detector sizes and therefore were renormalised accordingly.

5 Conclusions

We studied the correlations of particles in very small, extensive air showers.

We chose the method of scaled factorial moments, which seemed most appropriate for this type of analysis for a number of reasons:

  • Enables correlations to be determined over a wide range of the scale parameter \(\Delta \).

  • Makes it possible to focus on a particular type of correlation, building a set of “mixed events” that retains the effects that are not relevant in the study, but are in fact quantitatively dominant, which, when properly normalised, cancel them out.

  • Provides an opportunity to find specific patterns in the type of fractal structures that may or may not occur in cascading processes.

Calculations have shown strong non-Poissonian fluctuations of particle densities at small distances in the electromagnetic component of showers.

We have quantitatively determined the expected factorial cumulants on the basis of simulations with the CORSIKA program. In the range of primary energies studied, they are very strong for small distances. The correlations decrease rapidly with the size of the shower, the primary particle energy, as would be expected from the central limit theorem.

The “intermittent” behaviour of factorial cumulants (power-law dependence of \(F_k(\Delta )\)) is an interesting feature in itself and, rather surprisingly, is confirmed by experimental results.

Determining the value of the two-particle correlation is important for interpreting the results of shower measurements from very small, local shower arrays and measurements of the density spectrum of shower particles. The results obtained show that, for arrays of detectors located at distances of up to a few metres, the probability of two (three or even four) of them being hit simultaneously by a very small shower containing a total of several tens of particles is significantly higher than the trigger rate expected for a given shower particle density and a simple assumption that the shower particles reach the observation plane independently. The local density fluctuations are much larger than one would naively expect from a Poisson-like distribution. In small showers, where the particles are scattered over an area of a few or tens of thousands of square metres on average, the occurrence of a few particles at a small distance from each other, measured in metres, is more than 100 times more frequent than would be expected from a Poisson distribution with the same mean.

For large shower arrays designed to measure extensive air showers at the “knee” region and beyond, these effects are generally negligible, but small experiments and atypical systems sensitive to hits from single cosmic ray particles, where they can provide a messy background, should take these effects into account.