Feasibility of Correlated Extensive Air Shower Detection
with a Distributed Cosmic Ray Network

The global rate of cosmic ray shower events is quite substantial PDG (2018), and the chances of detecting physically-unrelated but nonetheless simultaneous showers increases with observatory size regardless of the technology employed. Two features of truly simultaneous showers which will help distinguish them from the combinatoric background are their opening angle $\Delta\psi$, which is expected to be small for nearly parallel correlated showers, and the great-circle separation between the two geographic locations, $\Delta s$, which tends to be larger for unrelated showers.

We assume that the incident flux of cosmic rays is isotropic, though some energy-dependent anisotropy in the UHECR flux has been identified with higher-order features Aab et al. (2015, 2014, 2017, 2018); Abbasi et al. (2014, 2018); di Matteo & Tinyakov (2018). Unrelated coincident showers are modeled by Monte Carlo methods, generating two random geographic locations drawn with uniform spherical density with random headings. The opening angle, $\Delta\psi$ is then computed from the difference in headings, and the great-circle separation between the two geographic locations, $\Delta s$, is computed from the well-known haversine distance formula.

The GZ effect is a potential mechanism for the generation of nearly-simultaneous, but greatly separated showers. UHECR nuclei are split by solar photons on their way to Earth, and the two daughter nuclei are separated by the solar magnetic field. Below, we describe the probability for photodisintegration, our model of propogation of the daughters in the solar magnetic field, and an efficient algorithm to generate simulated GZ dual showers. Knowledge of the heliospheric magnetic field (HMF) is the dominant source of systematic uncertainty in the entire analysis.

UHECRs striking the Earth’s atmosphere produce extensive air showers, which can be modeled with the the CORSIKA Heck et al. (1998) software package. As the computational time for a shower simulation scales roughly linearly with primary energy, running high-statistic simulations for ad hoc primary projectile, incident angle and observational altitude configurations demand significant computational costs. The statistical analysis of a search for rare events among a high-rate background demands a large number of sample showers. To ensure computational tractability, a parameterizated particle density as a function of the distance from the shower core is typically fit to samples of showers simulated with CORSIKA at various energies and for several species.

We compiled CORSIKA in 64-bit mode with GHEISHA 2002d for a horizontal flat array with thinning (including LPM) support. As the hadronic interaction model is the dominant uncertainty in Monte Carlo shower codes, we performed each simulation configuration 5 times, one for each available hadronic model: DPMJET-III (2017.1) with PHOJET 1.20.0, QGSJET 01C (enlarged commons), QGSJETII-04, SIBYLL 2.3c and VENUS 4.12. The charmed particle / tau lepton PYTHIA option was activated for DPMJET, QGSJET and SIBYLL. For all non-vertical simulations, the curved atmosphere selection was enabled. Values of essential CORSIKA parameters believed to optimize the realism of the shower simulations with acceptable computational time is provided in Tab. 2

For each of the 5 hadronic interaction models, 100 shower simulations were run for each of 5 shower-primary projectiles (photon, proton, Helium, Oxygen and Iron), for each of 8 shower-primary energies ($10^{14}$, $10^{15}$, $10^{16}$, $10^{17}$, $10^{18}$, $10^{19}$, $10^{20}$ and $10^{21}$ eV), and for each of 4 incident zenith angles (vertical, 30^∘, 60^∘ and 80^∘). For vertical simulations, 10 observation altitudes were specified at 0, 0.5, 1, 1.4, 2, 5, 10 and 20 km a.s.l. For angled-incident simulations where the curved atmosphere option limits observation altitudes to 1 per simulation, 5 simulations were performed for altitudes 0, 1, 2, 5 and 10 km (instead of all 5 hadronic models, only one was chosen at random per simulation). In total, roughly 100,000 simulations were performed.

While several standard paramaterization exist, it is not uncommon for both experiments and Monte Carlo simulations to find deviations from the analytically-motivated NKG Greisen (1956) distribution near and far from the core Barnhill et al. (2005); Hillas & Lapikens (1977); Fenyves et al. (1988). To accommodate these deviations, it is typical to develop variants of these parameterizations Gaisser & Hillas (1977); Schiel & Ralston (2007). Here, we develop a novel longitudinal shower density particle parameterization, which we name EKA, for photons and muons as a function of primary atomic mass, energy, incidence angle and observational altitude that describes well both the shower core and extended tails:

$$\rho(r;a_{n})=e^{a_{0}}r^{-a_{1}}exp\left(-\frac{r^{a_{3}}}{e^{a_{2}}}\right)$$

(7)

where the four $a_{n}$ coefficients are best-fit as a function of first-order factors of transformed primary mass number, $A$, primary energy, $\epsilon$, and altitude of observation, $h$ as:

$\displaystyle\begin{split}a_{n}(A^{*},\epsilon^{*},h^{*})\cong\;&c_{n}^{0}+\\ &c_{n}^{1}A^{*}+c_{n}^{2}\epsilon^{*}+c_{n}^{3}h^{*}+\\ &c_{n}^{4}A^{*}\epsilon^{*}+c_{n}^{5}A^{*}h^{*}+c_{n}^{6}\epsilon^{*}h^{*}+\\ &c_{n}^{7}A^{*}\epsilon^{*}h^{*}\end{split}$

(8)

where,

$\displaystyle\begin{split}A^{*}&=\mathrm{ln}\;(A+1)\\ \epsilon^{*}&=\mathrm{log_{10}}\;(\epsilon/10^{18})\\ h^{*}&=100\;\left(1-\mathrm{log_{10}}\;(10-h_{\mathrm{eff}}/10)\right)\end{split}$

(9)

with energy in electron-volts, altitude in kilometers and for vertical showers, $h_{\mathrm{eff}}$ is simply the altitude of observation, $h$. In general, for showers inclined by an angle $\theta_{0}$ above the observation horizon at altitude $h$, and coming from azimuthal direction $\phi_{0}$,

$$h_{\mathrm{eff}}=h^{\prime}+r\;\mathrm{sin}\;\theta_{0}\;\mathrm{cos}\;(\phi-\phi_{0})$$

(10)

where $h^{\prime}$ in turn is a function of the altitude of the first interaction, $h_{0}$:

$$h^{\prime}=h_{0}-(R_{E}+h)\left(\sqrt{\left(\frac{R_{E}+h_{0}}{R_{E}+h}\right)^{2}-\mathrm{cos}^{2}\;\theta_{0}}-\mathrm{sin}\;\theta_{0}\right)$$

(11)

where $R_{E}$ is the radius of the Earth in kilometers, and the average first interaction altitude, $h_{0}$, is modeled from simulation to good agreement as:

$$h_{0}(A^{*},\epsilon^{*})\cong d_{0}+d_{1}A^{*}+d_{2}\epsilon^{*}+d_{3}A^{*}\epsilon^{*}$$

(12)

with coefficients $d_{n}$ provided in Tab. 1. Example fits of our EKA model and the NKG model to CORSIKA samples are shown in Fig. 4. The EKA model is used to describe the expected particle densities in all subsequent simulations.

In this subsection, we describe our modeling of the effective observing power of a globally distributed array, and our methods for simulating dual showers impacting on the array, for two CRAYFIS array size and density scenarios.

The Monte Carlo model of the probability of randomly coincident cosmic rays is shown in Fig. 9 as a function of the opening angle, $\Delta\psi$, and the great-circle separation between their geographic locations, $\Delta s$. An analytical representation which agrees well with Monte Carlo results of $10^{7}$ simulated cosmic ray pairs is:

$$\begin{split}\mathrm{PDF}(\Delta s,\Delta\psi)=&\frac{1}{4R_{E}}\sin\left(\frac{\Delta s}{R_{E}}\right)\sin\left(\Delta\psi\right)\\ &\left[1+\frac{3}{4}\cos\left(\frac{\Delta s}{R_{E}}\right)\cos\left(\Delta\psi\right)\right].\end{split}$$

(15)

Constraints on the opening angle $\Delta\psi$ are the most powerful background discriminator as it is rare for random simultaneous showers to also be near-parallel. A model of the background such that $0<\Delta\psi<\Psi$, for some limit $\Psi$, Eq. 15 becomes,

$$\begin{split}\mathrm{PDF}(\Delta s;&\Delta\psi\leq\Psi)=\\ &\frac{1}{4R_{E}}\csc^{2}\left(\frac{\Psi}{2}\right)\sin\left(\frac{\Delta s}{R_{E}}\right)\\ &\left[1-\cos\Psi+\frac{3}{8}\sin^{2}\Psi\cos\left(\frac{\Delta s}{R_{E}}\right)\right]\end{split}$$

(16)

The methods described earlier are used to estimate the rate at which UHECRs photodisintegrate and produce daughters which both strike the Earth’s atmosphere.

As an example, the probability for an Oxygen UHECR nucleus to photodisintegrate is shown in Fig. 10 as a function of the distance from Earth. The average probability over all incoming trajectories to photodisintegrate and the subsequent fraction of these events that produce daughter particles which both strike the Earth is presented in Fig. 11a.

Photodisintegration via the Giant Dipole Resonance is most effective on nuclei with energy in excess of $10^{18}$ eV, generally plateauing somewhere between a 1 in 10 million likelihood for heavy elements, and 1 in 100 million for light elements. Somewhat conveniently, the gyroradius for fragments and protons at this energy range become large enough to almost ensure dual-EASs (Fig. 11b)—as the gyro-radius scales inversely with atomic number, $Z$; this is more the case for lighter elements than heavier.

The mean expected fraction of Earth-bound cosmic rays to produce dual-EASs (Fig. 12a) is the product of the fraction of UHECRs to photodisintegrate along their trajectory (Fig. 11a) and the probability to have both daugthers strike the Earth (Fig. 11b). The small probability for a GZ Effect candidate event is balanced by the very large size of the Earth and the high flux at lower energies, resulting in a large yearly global flux, as shown in Fig. 12b under the assumption that the cosmic rays are of a single elemental nucleus species.

As the GZ probability depends on the atomic number, a calculation of the overall rate of dual showers from GZ effects over the surface of the Earth must fold in the expected relative abundance of various nuclei. However, as this is not precisely known, each parent element’s flux contribution was estimated with relative weights that favor lighter elements (e.g. Oxygen, relative weighting $10^{-1}$) and suppress heavier elements (e.g. Uranium, relative weighting $10^{-3}$). Figure 12b shows that this weighting scheme favors GZ events from lighter element primaries with energies less than $10^{18}$ eV, and yields very slightly toward heavier primaries at the highest energies—a trend consistent with the observed cosmic ray spectrum PDG (2018). The observational consequence of this weighting scheme turning out to be opposite of reality would ultimately manifest itself in figure 18 as a left-translation (in minimum observation time) by an order of magnitude at most.

Additionally, the distribution of dual showers over the surface of the Earth is not isotropic, as it depends on the relative location of the Sun. A globally distributed sensor array which relies on consumer smartphones is likely to have a significantly higher effective area during the night, when phones are in darkness and not otherwise in use. Figure 13 shows the geographically binned expected flux per year for examples of higher and lower energy cosmic rays.

It can be seen that there is generally a day/night asymmetry—the highest-energy GZ parent UHECRs shown ($10^{20}$ eV) are nearly 4-times more likely to dual-EAS on the sunny-side of Earth than the dark-side. This asymmetry reverses, however, around $10^{17}$ eV where it becomes nearly 60-times more likely to dual-EAS on the dark-side of the Earth with $10^{15}$ eV UHECRs.

Unlike the flux distributions where the geographic bias changes with energy, the dark-side of the Earth consistently receives more closely-separated EASs than the sunny-side over all UHECR parent energies. The time-of-day asymmetries in flux and EAS separation do not come so much as a surprise considering the probability to photodisintegrate is greatest for “solar-passing” trajectories (Fig. 10) that experience highly Doppler-shifted (head-on) photons and high photon field density. Solar-passing trajectories also experience the greatest HMF strengths, which tend to cause greater separations. As the gyroradius (Eq. (5)) for individual products scales in proportion to energy (and inversely with atomic number), one or both low-energy products increasingly misses the Earth (especially on the sunny-side of the Earth) resulting in the expected flux inversion favoring the dark-side of the Earth at low energies.

On the other hand, although “Earth-direct” trajectories benefit from head-on (highly Doppler-shifted) photons, the photon field density is dramatically lower from this direction, resulting in photodisintegrations happening closer to the Earth on the average where the comparatively weak HMF (and shorter product propagation distance) results in more closely-separated EASs. The overall separation is greatest for the lowest energies, and smallest for the highest energy. An average over all geographic bins is provided in Fig. 14.

The expected rate of observed dual-EAS events in the CRAYFIS array under the two scenarios is the product of several factors: the rate of cosmic rays incident on the Earth, the probability for those particles to photodisintegrate and both decay products to strike the Earth, and the fraction of the resulting dual showers that are incident on CRAYFIS clusters of sufficient density. Given an UHECR particle which disintegrates, the average energy of the lower-energy decay product is $\approx 0.23E_{\textrm{UHECR}}$, given the energy splitting described in Eq. 6 under the imposed relative nuclei abundance in §III.2.

In addition, we consider selection requirements to reduce the rate of the combinatoric background. Correlated dual EASs from the disintegration of a single nucleus will arrive along parallel trajectories, while randomly coincident showers will be more broadly distributed in angle. Given an estimate of angular resolution of $\Delta\theta\approx 30^{\circ}$ and $\Delta\phi\approx 30^{\circ}$ from Ref. Whiteson et al. (2016), we estimate the efficiency of a ceiling on $\Delta\Psi$, see Fig. 15.

The background to signal rate drops, as expected, with increasingly tight ceilings on the opening angle. However, an overly restrictive selection will reduce sensitivity even as it enhances signal purity. To balance these competing effects, we require $\Delta\Psi<60^{\circ}$. The $\Delta\Psi$ selection efficiency is folded together with the cosmic ray flux so that the GZ probability and the efficiency to observe dual showers to estimate the dual-EAS flux in both scenarios, as well as the background rate, are shown in Fig. 16.

The expected rate of dual EASs from the GZ effect is quite low, due to the small probability of the GZ effect itself and the difficulty of observing dual showers. In scenario “P,” the rate of signal is many orders of magnitude below the observed flux. In scenario “U”, the large number of observing sites leads to an amplification of random coincidences, swamping the signal enhancement; however, this is the result of considering the entire world-wide network at once. Partitioning the globe into signal-optimized geographical regions of interest may show great improvement on background rejection. In fact, scenario “P” partially demonstrates this by restricting events to only those that approach from the dark side of the Earth, and the limited number of sensitive land patches implicitly apply a cut on possible shower separations. For scenario “U,” searches for events within fiducial bands of shower separation distances, primary energy, and geographic positioning relative the Sun are expected to provide powerful background suppressors; however this is left for future work.

We have used the GZ effect as a model for the expected signature of globally-coincident showers, but the larger question we address is whether a global array would be sensitive to any new, potentially unanticipated, source of simultaneous EAS. The power of a new observatory is to discover the unexpected; new physics scenarios could boost the rate of GZ disintegration, or alternative models could provide showers with similar global distributions. To address this broader question, we estimate the ability of the CRAYFIS array to observe coincident showers from some unknown process which creates GZ-like signatures at some higher rate.

The appropriate statistical task is then to estimate the fraction of observed coincident showers which come from a generalized rate-boosted GZ-like process, rather than random background. We have partially reduced the background rate by cutting on opening angle $\Delta\Psi$, but we have not yet taken advantage of the difference in expected distances $\Delta s$ between real and random coincidence. Figure  17 shows this is a powerful differentiator, favoring true coincidences at smaller separations due to the reduced combinations for random coincidences. The CRAYFIS sensitivity to GZ-like-effects that result in dual EASs is estimated using Monte Carlo pseudo-experiments with an unbinned likelihood in $\Delta s$.

The separation distances for each of $N_{\mathrm{tot}}$ simulated dual-EAS events are drawn from the combined signal and background models (Fig. 17), where $N_{\mathrm{tot}}$ is stepped by powers of 10 from $N_{\mathrm{tot}}=10^{1}$ to $10^{4}$, and for varying values of a rate-boosting factor, $B$, which simply increases the expected rate of showers relative to the GZ effect. The boost-factors for “Scenario U” range in powers of 10 from $B_{\mathrm{U}}=10^{0}$ to $10^{9}$, and for “Scenario P” from $B_{\mathrm{P}}=10^{0}$ to $10^{13}$—where the greater boost is capped to be consistent with the established UHECR flux. The minimal observation time to reject the background-only hypothesis at $3\sigma$ significance is shown in Fig. 18.