Open Questions in Cosmic-Ray Research at Ultrahigh Energies

1. Introduction

Cosmic rays with energy exceeding 10¹⁸ eV≡ 1 EeV, are referred to as ultrahigh-energy cosmic rays (UHECRs). Extensive air showers (EAS) produced when a UHECR interacts with an air nucleus in the upper atmosphere have been measured since their discovery by Pierre Auger in the 1930s. The first observation of an EAS with an energy of ~10²⁰ eV was made at Volcano Ranch in February 1962 (Linsley, 1963). The study of UHECRs has continued ever since, with increasingly large detector arrays. Nevertheless, many aspects of the nature of UHECRs remain an enigma: What is the origin of these particles? What is their mass composition? How do the astrophysical sources accelerate particles to such extreme energies?

This document summarizes the discussions that took place during the workshop “The High Energy Universe: Gamma-ray, Neutrino, and Cosmic-ray Astronomy” at the Munich Institute for Astro- and Particle Physics (MIAPP). We met for 1 month in March 2018 and had daily discussions and presentations about the status and future of the field of UHECR study. What have we learned about UHECRs in the last years? Which of the open questions can we expect to be able to address with forthcoming detector upgrades and proposed next-generation experiments? What are the requirements for probing remaining open questions and going forward in the study of UHECRs?

An overview of the current status of experimental measurements is given in section 2. Section 3 presents the open questions in the field of UHECRs. The theoretical models that successfully describe UHECR data are summarized. Predictions are given of the sensitivity of forthcoming and proposed experimental measurements to specific theoretical models and to the presented open questions in general. In section 4, upcoming and proposed Earth-based and space-based experiments are presented. We conclude in section 5, with our view of the outlook of the field, and a set of suggestions that we judge as beneficial for addressing open questions at ultrahigh energies in the coming years.

2. Status of Ultrahigh Energy Cosmic Ray Research

2.1. Anisotropy

The detection of an UHECR flux excess in the direction of a (few) prominent nearby source(s) would act as a pharos in the search for ultrahigh-energy accelerators. The volume of the Universe accessible at ultrahigh energies is limited by interactions with the extragalactic background light (EBL) and cosmic microwave background (CMB) to about 1 Gpc around 10¹⁹ eV, dropping down to a few hundreds of Mpc beyond 5 × 10¹⁹ eV (Greisen, 1966; Zatsepin and Kuzmin, 1966). As UHECRs are charged particles, their propagation is further affected by extragalactic and Galactic magnetic fields: the higher the rigidity (energy over charge), the smaller the deflection. Searches for UHECR anisotropies have consequently focused on large angular scales around 10¹⁹ eV, where the cumulative flux from multiple objects could possibly be seen despite magnetic deflections. At rigidities beyond ~ 20 EV, the trajectories of cosmic rays through extragalactic and Galactic magnetic fields are expected to be ballistic, with small (< 10 degree) deflections over 100 Mpc of propagation, motivating searches for small-scale anisotropies. Beyond this energy threshold, localized excesses at small (1°) to intermediate (30°) scales have been sought for (Aab et al., 2015d), possibly emerging from a few nearby objects.

Studies at large angular scales are often performed with ground-based observatories through a Rayleigh analysis (Linsley, 1975) in right ascension, α, of the UHECR arrival direction. Because of rotation of the Earth, the exposure of UHECR observatories only depends on declination, δ, when averaged over several years of observations. Using more than 8 years of full-operation data (12 years since the start of deployment), the Pierre Auger Collaboration discovered a modulation of the event rate in right ascension at E_Auger > 8 EeV with a post-trial significance of 5.4σ accounting for the search in two independent energy bins (Aab et al., 2017b).

Combining the right-ascension analysis with an azimuthal one, the anisotropy signal appears to be consistent with a dipolar modulation over ~85% of the sky covered by Auger. The amplitude of the dipole, $6.5_{-0.9}^{+1.3}\%$, is 10 times larger than that expected from proper motion in a cosmic-ray frame coincident with the CMB reference frame, suggesting an anisotropic distribution of UHECR sources within a few hundreds of Mpc. As shown in Figure 1, the direction of the dipole lies 125° from the Galactic center, disfavoring a Galactic origin for cosmic rays observed above eight EeV. This detection thus possibly constitutes the first observational piece of evidence for an extragalactic origin of cosmic rays beyond the ankle. Interestingly, further splitting events at E_Auger > 4 EeV into four energy bins, the Pierre Auger Collaboration found an indication at the 3.7σ level of growth of the dipolar amplitude with energy, expected from the shrinking horizon with increasing energy (Aab et al., 2018a). Given the sharp drop in statistics at the highest energies, searches for large-scale features remain under-constrained beyond E_Auger > 32EeV.

The Pierre Auger Collaboration has performed searches for intrinsic anisotropy at small angular scales at energies exceeding 40 EeV, by comparing the observed number of events within angular windows of a specified radius with that expected from an isotropic UHECR flux. The strongest excess revealed by this search is obtained at E_Auger > 54 EeV in a window of radius 12° centered on (α, δ) = (198°, −25°) (Aab et al., 2015d). This is shown on the left panel of Figure 2. Although the local significance obtained from this excess reaches 4.3σ, a penalization for the scan in energy and in search radius results in a post-trial value of 0.4σ (p = 69%). The Telescope Array (TA) Collaboration has performed a search for flux excesses at energies exceeding 10 EeV, 40 EeV and 57 EeV with 5 years of data. The largest excess, with a local significance of 5.1σ, is obtained at E_TA > 57 EeV in the direction (α, δ) = (147°, 43°) on a 20° angular scale (Tinyakov et al., 2016). The skymap of events with energy E_TA > 57 EeV, smeared on a 20° angular scale, is shown on the right panel of Figure 2. The colormap gives the pre-trials corrected significance of the observation in each direction in the sky with respect to background expectations. Accounting for the scan in search radius results in a penalized significance of 3.4σ, hinting at a possible over-density coined the TA “hotspot.” An update of the analysis presented with 7 and 10 years of data (Abbasi et al., 2018a; Sagawa, 2018) indicates no increase in the significance of the excess.

The directions with largest departures from UHECR isotropy have been compared with the position of nearby prominent objects. The two most significant excesses in the Northern and Southern hemispheres are located near the supergalactic plane, and multiple candidate sources have been discussed either within or outside from Collaborations. For example, in Fang et al. (2014) and Heet al. (2016), a ranking of gamma-ray emitting sources detected within 200 Mpc attempted to identify possible candidates for the TA hotspot, such as the starburst galaxy M82, blazars of BL Lac type such as Mrk 180 and Mrk 421, but also regular star-forming galaxies and galaxy clusters. Similarly, Cen A, an FR-I radio galaxy, or starburst galaxies such as NGC 4945 and M 83 have been pointed out as lying 10−20° away from the Southernmost significant UHECR excess. These sources are powerful X-ray and (or) γ-ray emitters and could potentially explain the UHECR flux from the TA hotspot region.

To reach a more complete view of the UHECR sky, cross-correlation studies against numerous astronomical catalogs have been performed within the Auger and TA collaborations, as well as by independent groups. Models often assume that the UHECR source distribution follows the distribution of luminous matter in the nearby Universe, based on radio — 3CRR catalog—or infrared—IRAS and 2MASS—or X-ray—Swift-BAT—or gamma ray— Fermi-LAT—observations. These models account for the expected energy losses and deflections of UHECRs during their extragalactic propagation (Kashti and Waxman, 2008; Oikonomou et al., 2013; Aab et al., 2015d; Tinyakov, 2018). While such studies have not yet revealed any statistically significant (> 5σ) departure from isotropy, a recent search against γ-ray bright sources, that accounted for their expected relative flux has unveiled an indication of excess UHECR flux at 4.0 σ post-trial in the direction of starburst galaxies (at E_Auger > 39 EeV), and at 2.7σ post-trial in the direction of jetted active galactic nuclei (AGN) at E_Auger > 60 EeV (Aab et al., 2018c). A search by the TA Collaboration with fixed parameters at E_TA > 43 EeV is consistent with the Auger result for starburst galaxies, but also with isotropy, indicating that the currently limited statistics from the Northern hemisphere is not sufficient to discriminate between the two hypotheses (Abbasi et al., 2018b).

2.2. Spectrum

Measuring the energy spectrum of UHECRs at high precision is of prime importance for understanding the origin and mechanisms of CR acceleration and propagation. Data at the highest energies have been accumulated for decades by AGASA (Yoshida et al., 1995), Yakutsk (Egorova et al., 2004), HiRes (Abbasi et al., 2008b), and more recently by the Pierre Auger Observatory and Telescope Array. Given the steeply falling energy spectrum, particularly above 5 · 10¹⁹ eV, event statistics is important. The statistical power of different observatories can best be compared by their integrated exposures. For illustration, after more than 20 years of operation, AGASA has reached an exposure of 0.18 · 10⁴ km² sr yr. As of ICRC 2017, the Telescope Array has collected 0.8 · 10⁴ km² sr yr, and Auger dominates with 9 · 10⁴ km² sr yr. This is a factor of 50 higher relative to AGASA and about a factor of 10 higher relative to TA and demonstrates the enormous progress that has been made during the last decade.

Both TA and Auger are hybrid observatories comprising a set of fluorescence telescopes and a surface detector array (Tokuno et al., 2011; Aab et al., 2015a). Their absolute energy calibration is based on the calibration of the telescopes and on knowing the fluorescence yield of the atmosphere. The details of the energy calibrations differ between the two observatories. Auger uses the absolute light yield and its wavelength dependence as measured by the Airfly Collaboration (Ave et al., 2008). TA uses the absolute yield measured by Kakimoto et al. (1996) at 337 nm and the wavelength dependence of the fluorescence yield measured by FLASH (Abbasi et al., 2008a). The dependence on atmospheric pressure, temperature, and humidity is treated identically by both collaborations using the reference model formula reported in UHECR2012 (Keilhauer et al., 2013). The corrections for the invisible energy of air showers are based on data in the case of Auger (Aab et al., 2019) and on Monte Carlo simulations in the case of TA. Finally, Auger uses the data-driven constant intensity method to account for the zenith angle dependence of shower absorption, while TA uses again Monte Carlo simulations (Ikeda, 2011). The joint working group of Auger and TA established that the relative differences between Auger and TA solely due these effects amount to 6 %. This is well in line with the total uncertainties of the absolute energy scales of 14 % in case of Auger (Verzi, 2013) and 21 % in case of TA (Abbasi et al., 2016). The additional contributions mostly stem from the absolute calibrations of the telescopes and from reconstruction methods (Dawson et al., 2013). Despite of these differences, a remarkable agreement in the energy scale of the two observatories is found up to about 10^19.4 eV. As demonstrated in Figure 3, (left) (The Pierre Auger and Telescope Array Collaborations, 2018), re-scaling the energy scale of each experiment by only 5.2 %, which is well within the aforementioned systematic uncertainties of the two experiments, provides an excellent agreement of their measured fluxes. However, above this energy, larger differences remain present, which cannot be accounted for by an independent scaling of their reconstructed energies. It will be important to understand whether this difference is caused by systematic uncertainties arising at the highest energies only, or whether it has an astrophysics origin related to seeing different parts of the sky. To study that question, the joint working group between the Auger and TA collaborations has generated energy spectra for the Southern sky, seen by Auger only, for the Northern sky, seen by TA only, and for the declination range −15° ≤ δ ≤ 24.8°, seen by both observatories. The energy spectrum for the common declination band is depicted in the right panel of Figure 3. Obviously, the agreement is much better, but some differences are still seen. It should also be noted that the energy spectrum measured by Auger does not show any significant declination dependence, but that of TA does. As it is still too early to draw definite conclusions about the source of the differences, the joint working group will continue their studies. It is also worthwhile to note that the declination dependence of the energy spectrum seen by TA should cause a significant anisotropy in the arrival directions of UHECR. This has been studied in Globus et al. (2017) and was found to be in tension with astrophysical models aimed at reproducing observational constraints on anisotropies.

Another important question related to the UHECR energy spectrum is about the origin of the flux suppression observed at the highest energies. The GZK cut-off was predicted 50 years ago independently by Greisen (1966); Zatsepin and Kuzmin (1966) and was claimed to be found by the HiRes collaboration in 2008 (Abbasi et al., 2008b). At the same time, the Auger collaboration reported a flux suppression at about the same energy and with a significance of more than 6σ (Abraham et al., 2008). Above 10^19.8 eV, TA has reported the observation of 26 events (Tsunesada et al., 2018) and Auger has reported 100 events (Fenu, 2017) by ICRC2017. However, these numbers cannot be compared directly due to the difference in the energy calibration of the experiments. We discuss more this problem in section 3.1.

2.3. Mass Composition

The most reliable technique to measure the mass composition of UHECRs is the simultaneous measurement of the depth, X_max, at which the number of particles in an air shower reaches its maximum and the energy, E, of the shower. These quantities can be directly observed with non-imaging Cherenkov detectors, radio arrays, and fluorescence telescopes. As of today, only fluorescence detectors have reached enough exposure to measure X_max at ultrahigh energies. After pioneering measurements from Fly's Eye (Baltrusaitis et al., 1985) and HiRes (Abu-Zayyad et al., 2000), the fluorescence technique is currently employed by the Pierre Auger Observatory (Abraham et al., 2010) and the Telescope Array (Tokuno et al., 2012). Traditional particle detector arrays are in principle also capable to estimate the energy and mass of cosmic rays, e.g., by measuring separately the number of muons and electrons at ground level, but usually with a worse resolution and, more importantly, larger theoretical uncertainties from hadronic interactions during the air shower development. The latter source of uncertainty can be eliminated by cross-calibrating the measurements with the X_max and energy of a subset of so-called hybrid events (air showers observed simultaneously with both, fluorescence and surface detectors).

The current data on the average shower maximum, 〈X_max〉, as a function of energy from fluorescence (Aab et al., 2014a; Abbasi et al., 2018d; Bellido, 2018) and surface detectors (Aab et al., 2017a) is shown in the left panel of Figure 4. The event-by-event fluctuations of the shower maximum, σ(X_max), are displayed on the right panel of Figure 4. Only the measurements with fluorescence detectors have enough resolution to determine the intrinsic (as opposed to detector-related) standard deviation of shower fluctuations. For comparison, the predictions of 〈X_max〉 of proton- and iron-initiated air showers simulations using hadronic interaction models (Ostapchenko, 2011; Pierog et al., 2015; Riehn et al., 2016) tuned to LHC data are shown as red and blue lines.

These measurements of the first two moments (mean and standard deviation) of the X_max distribution suggest that the composition of cosmic rays becomes lighter as the energy increases toward the ankle (until around 10^18.3 eV) and then becomes heavier again when approaching ultrahigh energies. The data points from the surface detector of Auger might indicate a flattening of this trend at ultrahigh energies, but more statistics are needed to confirm this finding. Note that, whereas 〈X_max〉 scales linearly with the average logarithmic mass of cosmic-ray primaries, a large value of σ(X_max) can either signify a light composition or a mixture of light and heavy nuclei, whereas a small value of σ(X_max) corresponds to intermediate or heavy composition with a small admixture of light elements (see, e.g., Linsley, 1983; Kampert and Unger, 2012).

For a more quantitive insight on the mass composition of UHECRs, the Pierre Auger Collaboration fitted templates of four mass groups (p, He, N, Fe) to the X_max distributions (Aab et al., 2014b; Bellido, 2018). The derived mass fractions are displayed in Figure 5 and reveal an interesting pattern of alternating dominance of certain mass groups. At low energies, there are hints for a rapidly disappearing contribution of iron, which is qualitatively in accordance with the “knee” in the flux of the heavy Galactic component at 10^16.9 eV reported by the KASCADE-Grande Collaboration (Apel et al., 2011). In addition to this heavy component, there seems to be a large fraction of intermediate-mass nuclei at low energy, possibly signifying a second Galactic component (Hillas, 2006; Thoudam et al., 2016). Above 10¹⁸ eV the flux of cosmic rays is dominated by light primaries. These have to be of extragalactic origin to avoid a large anisotropy toward the Galactic plane that would conflict with the level of isotropy of cosmic-ray arrival directions reported by Auger (Abreu et al., 2012a) and TA (Abbasi et al., 2017). As the energy increases, there is a trend that protons are gradually replaced by helium, helium by nitrogen, and there might be an iron contribution emerging above 10^19.4 eV when the statistics of the fluorescence measurement run out.Due to the limited statistics and unknowns about hadronic interaction models, this trend is still largely uncertain.

2.4. Neutral Secondaries: Ultrahigh Energy Photons and Neutrinos

Neutral secondaries including neutrinos and photons are expected to be produced when UHECRs interact with extragalactic background photons during intergalactic propagation. These secondary particles are also referred to as cosmogenic or GZK neutrinos and photons in the literature. Their flux mainly depends on the chemical composition, maximum energy of UHECRs, and the source evolution model (e.g., Takami et al., 2009; Kotera et al., 2010; see Hooper et al., 2005; Anchordoqui et al., 2008 for secondaries from heavy nuclei). In general, photopion production is more efficient than photo-disintegration in producing secondaries. Figure 6 presents the expected cosmogenic neutrino flux from Alves Batista et al. (2019), based on UHECR models that best fit the Auger spectrum and composition measurements.

Specifically, the orange shaded area covers the expectation of the best-fit models with 90% CL and assuming a source evolution following the AGN, star-formation rate (SFR), and γ-ray burst (GRB) redshift evolution (Alves Batista et al., 2019). The dark orange shaded area shows the neutrino flux of the best-fit scenario with 99% confidence level (CL). In this scenario, the source evolution is assumed to be a power law of the cosmic scale factor (1 + z)^m and the index m is left as a free parameter. Alves Batista et al. (2019) found that for energy spectral index between 1 and 2.2, their fit preferred negative source evolution, i.e., m < 0. This may be due to an actual evolution of sources, or an effect of cosmic variance and local over-density. In addition to Heinze et al. (2016), Romero-Wolf and Ave (2018), Das et al. (2018), Wittkowski and Kampert (2018), and Heinze et al. (2019) also predicted the cosmogenic neutrino flux based on fitting to the UHECR data.

Upper limits to the UHE neutrino flux have been obtained by the IceCube Observatory (Aartsen et al., 2018), the Auger Observatory (Bellido, 2018), and ANITA Gorham et al. (2018a). The blue shaded area shows the cosmogenic photons in the best-fit scenario of Alves Batista et al. (2019) and 99% CL. In more optimistic models, which assume larger maximum energy, 1–10 EeV photons may be observed. The gray shaded area presents such a flux, which covers the predictions by a range of models in Figure 7 of Decerprit and Allard (2011). For comparison, the upper limit of the differential photon flux in the bin of 10–30 EeV has been derived based in Aab et al. (2017c) and is shown as the blue solid line. For reference, we also present the high-energy neutrino flux measured by IceCube (Aartsen et al., 2016d; Kopper, 2018), cosmic rays (Abraham et al., 2008; Apel et al., 2013; Abbasi et al., 2018e), as well as the extragalactic gamma-ray background measured by Fermi-LAT (Ackermann et al., 2015).

The Figure 6, right shows the latest upper limits in searches for UHE photons. The strictest upper limits in this energy range come from Auger. The predicted cosmogenic photon fluxes are from Sarkar et al. (2011). With its current exposure, Auger constraints the photon fraction to be ≤ 0.1% above 10¹⁸ eV (Aab et al., 2015c; Aab et al., 2017c). Measurements with the Telescope Array surface detector provide complementary limits in the same energy range in the Northern Hemisphere (Abbasi et al., 2018c).

2.5. Hadronic Interactions at Ultrahigh Energies

Good understanding of hadronic multiparticle production is needed for being able to derive composition information from air-shower data. While measuring shower profiles using fluorescence and Cherenkov light allows an almost model-independent determination of the shower energy (up to a correction of the order of 10–15% for “invisible” channels Barbosa et al., 2004), there is no model-independent means for estimating the primary mass composition. The most productive approach is the detailed simulation of a library of reference air showers with Monte Carlo models that have been designed and tuned to describe hadronic multiparticle production at man-made accelerator experiments Engel et al. (2011). Hadronic interaction models of this type include EPOS (Liu et al., 2004; Werner et al., 2006; Werner and Pierog, 2007; Pierog et al., 2015), QGSJET (Kalmykov and Ostapchenko, 1989; Kalmykov et al., 1997; Ostapchenko, 2006a,b, 2014), Sibyll (Engel et al., 1992, 2017; Fletcher et al., 1994; Ahn et al., 2009; Riehn et al., 2016; Fedynitch et al., 2018), and DPMJET (Ranft, 1995; Roesler et al., 2001) for high-energy interactions, typically with a laboratory frame momentum larger than 100 GeV, and FLUKA (Ferrari et al., 2005; Böhlen et al., 2014) and UrQMD (Bleicher et al., 1999) for low-energy interactions. In general, a very good description of inclusive air-shower observables is obtained (see Abreu et al., 2011; Abu-Zayyad et al., 2012).

An important aspect of the hadronic interaction models is the extrapolation of accelerator data to center-of-mass energies of up to $\sqrt{s}~400$ TeV, well beyond energies accessible at colliders, to forward phase space regions not accessible in experiments, and to projectile and target particle combinations not measured in accelerator experiments. Given that we still cannot calculate predictions of QCD for the bulk of hadron production of importance for air showers, there is considerable ambiguity in modeling hadronic interactions. This ambiguity leads to model-dependent results for the mass composition as shown, for example, in Figure 5. Additional data from collider and fixed-target experiments and progress in the theory and phenomenology of multiparticle production are required to lower these uncertainties. For example, the LHC data at equivalent energy of $E_{\text{lab}}~10^{17}$ eV show a moderate rise of the proton-proton cross section and secondary particle multiplicity. Updating the interaction models led to a shift of the X_max predictions to larger depths (Engel et al., 2011; Pierog, 2018). While it was still possible to interpret the measured mean depth of shower maximum with a pure proton composition within the uncertainties using pre-LHC models, a mixed composition is clearly preferred if post-LHC models are applied. The shower-by-shower fluctuations of X_max provide an even stronger constraint; see Figure 4. The depth of the first interaction point of an air shower is exponentially distributed, dP/dX₀ ~ exp(−X/λ), with λ the interaction length. Hence, the fluctuations of X₀ are σ(X₀) = λ. Using the measured values of the proton-air cross section (see Figure 7, left), one gets σ(X₀) ~ 50 g/cm². Even if there were no additional fluctuations introduced by the shower evolution from the first interaction to the shower maximum, the proton-air cross section would have to be two times larger to bring these fluctuations down to 25 g/cm². Such a drastic increase in the proton-air cross section would violate unitarity constraints in QCD and would require a new type of interaction taking over at energies beyond 2 × 10¹⁸ eV.

Air-shower measurements can also be used to derive information on hadronic interactions. Given that the primary cosmic-ray composition appears to be mixed in the energy range of relevance here, there is typically a strong correlation between the results of such measurements and the assumed primary mass composition. An exception is the measurement of the proton-air cross section. If done in an energy range in which there is a large fraction of protons in the mass composition of cosmic rays, one can select showers that develop very deep in the atmosphere to build a proton-dominated sample. Then the depth fluctuations can be related to the proton-air cross section for particle production. Recent results are shown in Figure 7, left.

There is increasing evidence for a discrepancy between the number of muons predicted by model calculations and that measured at very high energy. One of the most direct measurements demonstrating this muon discrepancy is shown in Figure 7, left. Depending on the interaction model used for reference and the measurement, there are about 30 − 60% more muons found in data than predicted. This muon puzzle is one of the most important problems in hadronic interaction physics as it is very difficult, if not impossible, to increase the number of muons by such a large fraction just by changing the physics of the first interaction. Enhanced production rates of baryon-antibaryon pairs (Pierog and Werner, 2008) and ρ⁰ mesons in air showers (Drescher, 2008; Ostapchenko, 2013) have been shown to have a large impact on the muon number. While NA61 measurements (Aduszkiewicz et al., 2017) have confirmed an enhanced forward production rate of ρ⁰ mesons, no increased proton-antiproton production rate has been found at LHC. Even though tuning these production processes increases the predicted muon number (Riehn et al., 2016), the discrepancy to air shower measurements still persists. It is likely that not only the number but also the production depth (Aab et al., 2014c; Collica, 2016), energy spectrum and, hence, the lateral distribution of muons is not well described by the models.

3. Open Questions

3.1. Precision Measurements of Spectrum and Mass-Composition

3.1.1. Relevance of the Energy Resolution

Enormous progress has been made recently from observing simple all-particle power-law distributions with just seeing the knee and ankle of the cosmic-ray spectrum, to uncovering a much more complex structure with an additional “second knee” at about 10¹⁷ eV, an ankle-like structure between the knee and this second knee, and the steep cut-off at the highest energies. Moreover, not only all-particle spectra can be derived from the air-shower data, but also energy spectra of different mass groups. All these achievements provided new insight into the astrophysics causing those structures. This became possible only by advancing both the precision of air-shower observations and reconstructions and the statistics of the data. In fact, improving simultaneously the quality and quantity will also be the key to making progress in the future.

The disentanglement of the all-particle energy spectrum into that of individual mass groups from about 10¹⁵ eV to 10¹⁷ eV, most notably by KASCADE and KASCADE-Grande, has provided new insights into the origin of the knee and ankle and will be discussed in section 3.2.2 in the context of the transition from Galactic to extragalactic cosmic rays. The origin of the flux suppression of cosmic rays at highest energies is still debated. The two competing explanations are energy-losses of UHECR in the CMB or nearby sources of UHECR with corresponding maximum acceleration energies (see section 3.2.5).

Identifying individual sources of UHECR would answer that question and remains the ultimate goal of future studies. It can be expected that the arrival directions of light primaries at the highest energies are correlated with UHECR sources located within the GZK sphere. Identifying such sources calls for a precise shower-by-shower determination of the energy and mass of the primary particle to avoid cosmic rays of lower energy diluting the event sample and to avoid heavy primaries, suffering stronger deflections, blurring the source spots in the sky. In both cases, the experimental energy and mass resolution determine the ratio of possibly source-correlated events to background events so that compromises in experimental resolution need to be paid for by larger event statistics, i.e., by larger exposures. Obviously, the steeper the spectrum in the region of interest, the stronger is the effect of spillover. This has been studied in a simplified model in Brümmel et al. (2013), depicted in Figure 8. Here, the blue line represents the energy spectrum observed with 30 % energy resolution from the true parent distribution (shown as dashed line). The red line shows the distribution of events that leak into the region of interest (above the vertical black line) despite having a true energy $E_{\text{true}}<E_{\text{thresh}}=10^{19.63}$ eV. Assuming all events below E_thresh being isotropic and those above being correlated to sources, 58 % of the events observed above the applied threshold would be isotropic background and dilute the signal. To compensate for this unwanted effect, the applied energy threshold could be increased as is illustrated in the right panel of Figure 8. Increasing the threshold of a 30 % energy resolution detector from 10^19.63 eV to 10^19.83 eV would yield the same “signal purity” as a 10 % energy resolution detector has reached at 10^19.63 eV. However, the flux of UHECRs drops at the same time by an increase of a factor of two, so that the relaxation in energy resolution from 10 % to 30 % needs to be compensated by a factor two in exposure. We understand this to be a simplified model, but it serves the general discussion and demonstrates the importance of the effect.

3.1.2. Composition at Ultrahigh Energies

As mentioned in section 2.3, the inferred cosmic-ray composition at Earth shows a peculiar dependence on energy (cf. Figure 5). The sequence of alternating groups of elements and the increase of mass with energy could be caused by a Peters cycle (Peters, 1961) at the accelerators, i.e., maximum energy that depends on rigidity R = E/Z or due to photonuclear spallation processes during propagation (e.g., Allard et al., 2007, 2008; Hooper and Taylor, 2010) to Earth and in the source (Unger et al., 2015), leading to scaling with energy per nucleon E/A.

The factor between the maximum fraction of protons and helium in Figure 5 is close to 4, which would favor a spallation scenario. However, the data does not yet constrain the maximum of the N and Fe group and, moreover, a combination of Peters cycle and spallation effects is not excluded. Only the detection of several cycles (if any) will allow for an unambiguous disentanglement of the combined effect of spallation effects from the propagation to Earth and the possible existence of a Peters cycle and/or photonuclear interactions in the source. For this purpose, large-exposure observatories with a good (equivalent or better to current fluorescence detectors) mass resolution is needed.

Another important open question related to mass composition is the evolution of the rigidity R = E/Z with energy. The angular deflections are proportional to θ ∝ 1/R and for an ensemble of different charge-groups with fraction f_i and charge Z_i it is 〈θ〉 ~ 〈1/R〉 = ∑f_iZ_i/E. The evolution of the average rigidity with energy is shown in Figure 9. As can be seen, the rigidity is increasing with energy and therefore the angular deflection in magnetic fields should decrease with energy, i.e., the increase in the average mass of cosmic rays with energy as shown in Figure 5 is slow enough to not outrun the increase of energy. No high-quality data currently exists at ultrahigh energies where hints for anisotropies at intermediate scales were reported. Note that the average logarithmic mass derived from 〈X_max〉 is not enough to determine the rigidity, because the mass-to-charge ratio is 1 for protons and ~ 2 for other elements.

3.2. Astrophysics

3.2.1. Origin of the Bulk of UHECRs

The challenge of accelerating cosmic rays to 10²⁰ eV was succinctly presented in the form of the minimum requirement for the accelerators, in what is now commonly referred to as the “Hillas condition” (Hillas, 1984). It states that a necessary condition to accelerate particles to ultrahigh energy is that of confinement; particles can stay in the acceleration region as long as their Larmor radius is smaller than the size of the accelerator. Thus, the maximum energy achievable, E_max, in a source with characteristic size, R, and magnetic field strength, B, is, E_max = eBR. Here, R = l · Γ, with l the comoving size of the source, and Γ the Lorentz factor of the motion, which is thought to be Γ ~ 10 − 50 in AGN jets (e.g., Lister et al., 2019), and Γ ~ 10 − 1, 000 in GRBs.

For acceleration in a shock with velocity β_sh (in units of the speed of light), the maximum achievable energy is,

where η parametrises the efficiency of acceleration, with η = 1 the maximum achievable efficiency when diffusion proceeds in the Bohm limit.

The confinement condition is not sufficient to guarantee cosmic-ray acceleration to 10²⁰ eV. This depends on the details of the acceleration mechanism and the timescale for energy loss in the source environment. A summary of constraints on astrophysical sources based on the Hillas condition was presented in Ptitsyna and Troitsky (2010).

Figure 10 shows classes of objects in terms of the product of their radial size, R, magnetic field strength, B, and associated uncertainty in the ideal limit where η = 1. The solid diagonal lines show the minimum product of BR required to accelerate protons (red) or iron nuclei (blue) to 10²⁰ eV for a fast shock where β_sh = 1. Classes of objects to the left of the lines do not satisfy the Hillas criterion. As shown with the dashed diagonal lines, the required product of BR is higher for slower shocks (β_sh = 0.01 is shown for illustration). The plot reveals that normal galaxies, supernovae, and stars that drive massive magnetized winds such as Wolf-Rayet stars do not satisfy the confinement condition. For the other source classes in the plot, the confinement condition is satisfied.

Another condition that must be met by UHECR accelerators is that they must possess the required energy budget to produce the observed UHECR diffuse flux. The energy production rate of UHECRs has been estimated in Waxman (1995b), Berezinsky et al. (2006), Katz et al. (2009), and Murase and Takami (2009) under the assumption that UHECRs are extragalactic protons. Most recently the energy production rate of UHECRs was estimated in Aab et al. (2017d), where a combined fit to the all-particle spectrum and X_max distributions at energy 5 × 10¹⁸ eV and beyond measured at the Pierre Auger Observatory was performed. Here, a mixed injected composition was allowed. The best-fit model corresponds to a UHECR energy-production rate, $E_{\text{UHE}}{Q}_{E_{\text{UHE}}}\approx 5\times 10^{44}\text{ erg }{\text{s}}^{-1}{\text{yr}}^{-1}$. The true value of the UHECR energy budget depends on the source-by-source injected spectrum, composition, and luminosity density evolution of the sources, and may differ from that of Aab et al. (2017d). Further, the inferred UHECR production rate depends on the chosen energy range (see e.g., Murase and Fukugita, 2018). Most estimates converge to $E_{\text{UHE}}{Q}_{E_{\text{UHE}}}\gtrsim 5\times 10^{43}\text{ erg }{\text{s}}^{-1}{\text{yr}}^{-1}$.

Figure 11 shows the energy budget of various source classes based on infrared, radio, X-ray, and gamma-ray observations, and compares it to the UHECR production rate estimated in Aab et al. (2017d). We used characteristic luminosities for each source type and the luminosity density at z = 0, motivated by the fact that locally observed UHECRs must originate in nearby sources located at ≲ 100 Mpc. The solid diagonal line shows the required energy budget to power observed UHECRs assuming that the UHECR luminosity of the sources, L_cr, is equal to the luminosity of the sources in the wavelength studied, L_γ. Sources to the left of the line do not satisfy the energy budget condition. The UHECR luminosity of individual sources need not be equal to their radiative luminosity. In the absence of knowledge of the relation between the two, we show for illustration dashed diagonal lines for the condition L_CR = 0.1 × L_γ and L_CR = 10.0 × L_γ. Note that the Hillas criterion imposes an independent lower limit on the magnetic luminosity of a UHECR source as shown in Equation (4) [see relevant discussion in Lemoine and Waxman (2009); Fang and Kotera (2016)].

The orange dashed line gives the minimum source number density constraint, which comes from the analysis of arrival directions of UHECRs detected in Auger of Abreu et al. (2013). The lack of significant clustering in the arrival directions of UHECRs with energy exceeding 70 EeV was used to derive a lower limit to the UHECR source number density, considering that UHECRs might have suffered deflections as large as 30°. Classes of steady sources to the left of the orange diagonal line do not satisfy the source number density constraint, unless UHECR deflections are significantly larger than investigated in Abreu et al. (2013).

In order to compare the energy budget constraint to the energy budget of transient source classes, the observed burst rate, ρ, must be converted to the effective number density for UHECRs, n_eff = (3/5)ρ · τ (see e.g., Murase and Takami, 2009), with τ the apparent burst duration of the UHECR burst (Waxman and Miralda-Escude, 1996),

where D is the distance traveled by the UHECR, λ the correlation length of the regular magnetic field, and Z the atomic number of the UHECR nucleus. Similarly, the effective luminosity can be estimated by modulating the burst fluence by τ. Figure 11, also shows the effective luminosity and number density for transient sources, where we have used τ = 3 × 10⁵ yr. The effective number density shown for bursting sources is only valid for the assumed value of τ, which corresponds to mean extragalactic magnetic field strength 1 nG. Stronger magnetic fields would imply larger τ and hence, larger effective number density. In this case, the bursting sources satisfy the number density constraint more comfortably, but the effective luminosity also decreases so the comparison with the energy-budget constraint does not change. On the other hand, τ cannot be arbitrarily small. A lower limit comes from the time spread induced from the coherent component of the Galactic magnetic field, τ_min ~ 300 − 3, 000 years (see Murase and Takami, 2009 for details).

Below, we discuss the most plausible UHECR-source candidates in turn.

Gamma-ray bursts and energetic supernovae. Gamma-ray bursts are during their short lives some of the most spectacularly bright objects in the sky. They have long been discussed as likely sites of UHECR acceleration (Vietri, 1995; Waxman, 1995a). In general, GRBs are thought to easily satisfy the maximum energy requirement (see however Samuelsson et al., 2018). Inspection of the energy budget diagram reveals that high-luminosity GRBs are roughly consistent with the energy budget requirement, though on the low side. As cautioned earlier, the UHECR energy budget is uncertain and consistent with being ten times lower than the model shown in Figure 11 as a benchmark.

Low-luminosity GRBs, which are a less-well-known source population, seem to occur with a much larger rate locally than high-luminosity GRBs. They are appealing as sources of UHECRs (Murase et al., 2006; Liu et al., 2011) as the relatively milder radiation fields with respect to those of high-luminosity GRBs would better allow the survival of UHECR nuclei. Several articles have addressed the conditions of acceleration and survival of nuclei in high-luminosity (Murase et al., 2008b; Wang et al., 2008; Horiuchi et al., 2012; Baerwald et al., 2015; Globus et al., 2015), and low-luminosity GRBs (Murase et al., 2008b; Horiuchi et al., 2012; Boncioli et al., 2018; Zhang et al., 2018) and find regimes in which GRBs could power all the observed UHECRs and be consistent with the UHECR composition measurements.

Though standard supernovae are not expected to be able to accelerate cosmic rays to ultrahigh energies, the ejecta of trans-relativistic and engine-driven supernovae which typically reach mildly relativistic speeds may also be able to accelerate UHECRs (Wang et al., 2007; Chakraborty et al., 2011; Liu and Wang, 2012; Fang et al., 2018; Zhang and Murase, 2018). A feature of GRB and engine-drive SN models is that the composition resulting from stellar evolution models can explain the UHECR composition data observed by Auger.

In 2017 the detection of gravitational waves from the merger of a neutron star binary, followed by a short GRB and electromagnetic emission from the remnant marked the discovery of this, long-sought-for, class of events (Abbott et al., 2017a,b). In Rodrigues et al. (2019) it was shown this class of sources could be producing the cosmic rays observed right below the ankle. On the other hand, Kimura et al. (2018a) showed that the tail of the Galactic cosmic-ray spectrum around the second knee can be explained by remnants of Galactic neutron star mergers.

A brief mention to the winds of Wolf-Rayet stars is also due here. Though inspection of Figure 10 reveals that the winds of these sources likely do not satisfy the Hillas criterion for 10²⁰ eV UHECRs, (see however Biermann and Cassinelli (1993) for a different view) the magnetized, powerful winds they drive have been proposed as possible acceleration sites of cosmic rays up to 10¹⁸ eV and could thus be responsible for the end of the Galactic cosmic-ray spectrum (Thoudam et al., 2016; Murase and Fukugita, 2018).

Active galactic nuclei. Active galactic nuclei (AGN) with powerful jets have long been considered as promising candidate sources of UHECRs. AGN with jets pointing to the Earth, referred to as blazars, would be the natural candidates if UHECRs escape the sources beamed and do not suffer severe deflections (Dermer and Razzaque, 2010; Resconi et al., 2017; Rodrigues et al., 2018). The signature of UHECR acceleration could be detectable in the gamma-ray spectra of blazars (Gabici and Aharonian, 2005; Essey and Kusenko, 2010; Essey et al., 2010; Kotera et al., 2011; Murase et al., 2012; Prosekin et al., 2012; Aharonian et al., 2013; Takami et al., 2013; Oikonomou et al., 2014; Tavecchio, 2014). However, the present-day density of nearby blazars shown in Figure 10 suggests that blazars alone do not satisfy the number density constraint. On the other hand, radio galaxies, the parent population of blazars (BL Lacs and FSRQs) with jets pointing away from the line of sight, are also UHECR source candidates, with Cen A, the nearest radio galaxy, a long standing candidate (Rachen and Biermann, 1993; Romero et al., 1996; Atoyan and Dermer, 2008; Dermer et al., 2009; Gopal-Krishna et al., 2010; Biermann and de Souza, 2012; Wykes et al., 2018). In recent literature, several models have been proposed, which show that the observed UHECR flux and composition can be produced by radio-galaxies under different assumptions about the acceleration mechanism at the sources, namely shear (Kimura et al., 2018b) and “one-shot” re-acceleration (Caprioli, 2015). The re-acceleration models can explain the nucleus-rich composition data observed by Auger.

In jetted AGN, a lot of the power goes to energizing the lobes, which are very extended features with relatively small magnetic fields (B ~ 10⁻⁵ G) and proposed sites of UHECR acceleration (Takahara, 1990; Rachen and Biermann, 1993). It was recently shown in Matthews et al. (2019) with hydrodynamical simulations that acceleration to 10²⁰ eV is possible in these regions, and in Eichmann et al. (2018); Matthews et al. (2018) that nearby radio galaxies are strong UHECR candidates.

In addition, radio-quiet, low-luminosity AGN and quasar outflows have been discussed as possible sources of UHECRs (Pe'er et al., 2009; Duan and Caramete, 2015; Wang and Loeb, 2017). These are less powerful individually than jetted AGN but significantly more numerous.

Tidal disruption events. Stars that pass within the tidal radius of a super-massive black hole are disrupted and a large fraction of the resulting debris gets accreted onto the black hole. If the disruption occurs outside the black hole horizon a luminous flare of thermal emission is emitted and in a fraction of these events a jet forms (Hills, 1975; Rees, 1988). Only a handful of jetted TDEs have been observed to date, whereas the total number of known and candidate TDEs is at present close to 100. It was shown in Farrar and Piran (2014), based on the analysis of the prototypical jetted-TDE Swift J144+57, that jetted-TDEs can likely produce the bulk of observed UHECRs. The expected UHECR output from TDEs was more recently studied in Zhang et al. (2017), Dai and Fang (2017), Biehl et al. (2018), and Guépin et al. (2018) in the internal shock model. The above analyses conclude that given the relatively low inferred rate of jetted TDEs based on Swift data, whether the energy-budget constraint is satisfied depends intricately on the relation between the TDE radiative luminosity and UHECR luminosity. Based on theoretical arguments Farrar and Piran (2014), Biehl et al. (2018) showed that the energy-budget constraint is likely satisfied, despite the apparent failure of TDEs to satisfy the constraint based on the Swift data as shown in Figure 11.

Intermediate-mass black holes may also tidally disrupt stars. Depending on the combination of masses of both objects, tidal squeezing may trigger nuclear burning in the core of white dwarfs, leading to a supernova and potentially accelerating cosmic rays to ultrahigh energies (Alves Batista and Silk, 2017; Zhang et al., 2017; Guépin et al., 2018).

Starburst galaxies. Starburst galaxies are galaxies that are undergoing intense star-formation activity, typically demonstrated by infrared luminosities > 10 times higher than normal galaxies. They are observed to drive powerful, magnetized “winds” (nuclear outflows), which might be sites of high-energy particle acceleration (Anchordoqui et al., 1999). The maximum UHECR energy that can be achieved in the wind driven by starburst galaxies was recently studied in Anchordoqui (2018); Romero et al. (2018) and Murase and Fukugita (2018), with conflicting conclusions as to the feasibility of UHECR acceleration in starburst winds. Another natural possibility is that UHECR acceleration can occur in the disproportionately frequent extreme explosions that take place in starburst galaxies due to the high star-formation activity. These include low-luminosity gamma-ray bursts, trans-relativistic supernovae, and hypernovae, which do not have to occur only in low-metallicity environments (Zhang and Murase, 2018; Zhang et al., 2018).

Galaxy clusters. Galaxy clusters, the largest bound objects in the Universe, have also been considered as possible sites of UHECR production (Kang et al., 1997; Ryu et al., 2003; Murase et al., 2008a; Kotera et al., 2009). Though they possess moderate magnetic fields ~ μG (see e.g., Kim et al., 1991; Carilli and Taylor, 2002) they are extremely extended ~ 2 − 3 Mpc, and should thus be able to confine particles to extremely high energies (Fang and Olinto, 2016). Galaxy clusters could otherwise act as “reservoirs” which contain sites of UHECR acceleration, for example, jetted AGN (Murase et al., 2008a; Kotera et al., 2009; Fang and Murase, 2018).

Pulsars. Pulsars, the smallest and most highly-magnetized objects shown in Figure 10, induce strong magnetic potentials that can potentially also accelerate UHECRs (Gunn and Ostriker, 1969; Blasi et al., 2000; Arons, 2003; Murase et al., 2009; Philippov and Spitkovsky, 2018). Since they are the product of the death of massive stars and shrouded by a remnant enriched in heavy elements, it has been shown that they may produce UHECRs rich in nuclei (Fang et al., 2012, 2013).

3.2.2. Galactic to Extragalactic Transition

The cosmic-ray spectrum features three distinct spectral breaks in the energy range between 10¹⁵ and 10¹⁸eV. In order of increasing energy, these are the “knee,” “second-knee” (or “iron-knee”) and “ankle,” illustrated in Figure 12. Below we discuss the origin of each of the three features and viable scenarios for the transition between Galactic and extragalactic cosmic rays.

The physical origin of the knee feature remains unclear. Both a propagation and sources maximum energy origin of this feature remain viable candidates. An outline of these two scenarios is given below.

Propagation origin of the knee. Galactic cosmic rays are believed to diffuse within the Galactic magnetic turbulent sea. Within the plane of the Galactic disk, the dominant drivers of this MHD turbulence are believed to be supernova remnants (SNR), which inflate bubbles tens of parsec in size, driving magnetic turbulence on this scale (λ_max). Although turbulence is driven on such a scale, it subsequently cascades down to smaller wave modes, eventually terminating at the dissipation scale, λ_min.

Assuming that cosmic rays of a given Larmor radius r_Lar predominantly scatter resonantly from magnetohydrodynamic (MHD) turbulent modes λ of the same size (i.e., r_Lar = λ), the multi-PeV energy scale denotes the energy range at which the abundance of such modes diminishes rapidly. PeV cosmic-ray protons in ~ μG Galactic magnetic fields possess Larmor radii of 1 pc. The PeV energy scale is therefore motivated to denote the energy range in which CR diffusion within the Galactic disk magnetic field becomes inefficient (i.e., r_Lar ~ λ_max) (Giacinti et al., 2014). In such a scenario, the confinement of cosmic rays at higher energies becomes significantly less efficient, giving rise to a steepening of the cosmic ray spectrum capable of explaining the shape of the knee feature. More generally, any scenario in which a change in the transport regime occurs leading to inefficient confinement can also lead to this feature ( e.g., the transition from diffusive propagation to particle drift Ptuskin et al., 1993).

It is important to note, however, that such propagation origin scenarios for the knee make a considerable implicit assumption. For these scenarios it is necessary that luminous Galactic cosmic ray sources exist, capable of accelerating particles to energies well beyond the knee energy. Within the framework of our current understanding of Galactic SNR accelerators, however, such an assumption presents a considerable challenge (Bell et al., 2013). Indeed, presently, the only known Galactic source capable of achieving acceleration to the PeV scale is Galactic nucleus, Sgr A* (Abramowski et al., 2016), whose cosmic ray luminosity at these energies appears to be rather low (see e.g., Fujita et al., 2017).

Maximum-energy origin of the knee. Alternative to this propagation origin of the knee is the possibility that the PeV energy scale denotes the maximum energy of their Galactic sources, believed to be SNRs. An application of the Hillas criterion in Equation (1) to SNR gives a maximum energy of,

where the factor η describes how close to Bohm diffusion the maximum energy particles in the source achieve, β_sh is the shock velocity in units of c, B is the magnetic field in the acceleration region and R is the size of the source. Equation (3) indicates the need for considerable magnetic field enhancement to occur in order for such sources to act as effective PeVatron candidates. Such an enhancement may occur by the Bell mechanism (Bell, 1978a,b) in which CR accelerated by the SNR run ahead of the shock, whose current drives an instability in the upstream medium enhancing the upstream magnetic field present. Furthermore, observationally, there is now growing evidence that such magnetic field enhancement takes place within these sources. However, whether SNRs are actually able to accelerate up to the knee energy (3 PeV) remains an open question.

In either the propagation or maximum-energy scenario which describes the origin of the knee feature at 3 PeV, a family of corresponding knee features for the other nuclear species are naturally expected. Observationally, it remains unclear whether the composition of CRs at the energy of the knee feature (3 PeV) are protons, helium, or heavier species. Assuming the composition of the knee to be dominated by protons (i.e., a proton knee at 3 PeV), a corresponding iron knee feature at 100 PeV would be expected. Observational evidence for such a second knee feature was reported from the analysis of the KASCADE-Grande data (Apel et al., 2011).

On theoretical grounds, it remains extremely challenging for known Galactic CR accelerators to accelerate protons above PeV energies. The known magnetic field amplification scenarios place a hard cap for maximum energies achievable by SNR (Bell et al., 2013).

In addition, the low level of anisotropy of cosmic rays in the energy range 10¹⁷ − 10¹⁸ eV also disfavors a Galactic origin of any light component in this range (Abreu et al., 2012; Giacinti et al., 2012).

At energies at/just above that of the second knee, observational evidence suggesting the onset of a new component in the light composition spectrum is found in the KASCADE-Grande data, referred to as the proton ankle. Evidence pointing in this direction is also supported by the low-energy Auger HEAT X_max data, which show a lightening in the composition above 100 PeV. If the interpretation of these observational results is correct, the onset of this new light component marks the beginning of an extragalactic component in the arriving CR flux. Such an interpretation has considerable implications, which provide the possibility to shed new light on the extragalactic origins of these protons.

Extragalactic cosmic-ray protons at EeV energies undergo frequent Bethe-Heitler energy loss interactions with CMB photons, losing their energy through this process on Gyr timescales. These losses give rise to electron/positron pairs, $p+{\gamma }_{\text{CMB}}\to p+e^{+}{e}^{-}$, which subsequently feed electromagnetic cascades, with the energy flux cascading down to energies below 100 GeV, contributing to the diffuse gamma-ray background, particularly since such a source evolution allows for a Fermi type source injection spectrum (Taylor et al., 2015). Recent improvements in our understanding of the contributions to this background constrain the allowed level of these losses, which could prefer scenarios with negative evolution, i.e., that these extragalactic cosmic rays have a small filling factor in extragalactic space (Liu et al., 2016) (but see also e.g., Fang and Murase, 2018).

3.2.3. Source Identification Beyond the Ankle

The large-scale anisotropy discovered beyond the ankle by the Pierre Auger Observatory appears to be consistent with the distribution of extragalactic matter traced by near-infrared observations from 2MASS (Aab et al., 2017b). This can be seen as the first observational evidence, also supported by theoretical expectations, for UHECRs beyond the ankle originating from extragalactic sources. Most likely, not all galaxies behave as UHE accelerators, so that the question of which galaxies or galaxy types host UHE accelerators remains open.

Cross-correlation with catalogs of objects observed throughout the electromagnetic bands has proven a powerful means to address the question of possible associations. Such searches recently hinted (3–4σ) at a fraction of 10–15% of UHECR events being consistent with the directional and flux distributions expected from either extragalactic matter—traced by 2MASS or Swift-BAT X-ray observations—or specific types of extragalactic sources—starburst galaxies and jetted AGNs—traced by their radio and gamma-ray emission (Abbasi et al., 2018b). Even if such an anisotropic signal reached the 5σ discovery threshold in the near future, it probably would not be sufficient to claim identification of UHECR sources. As correlation does not imply causation, a necessary condition for an identification of some or all the sources would be to leave as little room as possible for a confounding variable, that is, a hidden variable causing a spurious correlation. Such a feat would require coverage over the entire celestial sphere—to avoid blind regions where a different source type could contribute—constraints on the redshift evolution of the UHECR production rate—to enable a tomographic probe of source populations—and completeness in the source models up to the propagation horizon and down to a sufficiently low luminosity.

Ground-based observations come with a partial view of the celestial sphere. Nonetheless, attempts at full-sky coverage by combining data from the largest Northern and Southern observatories have been performed by the Pierre Auger and Telescope Array Collaborations. Such an approach is limited by the mismatch in energy scale between the two experiments, which could cause spurious anisotropies due to an improper contrast between the flux inferred from each dataset. The collaborations have designed a method to match the flux in the declination band covered from both sites, providing a common view on the UHECR sky beyond the ankle and above the flux suppression (Aab et al., 2014d; Biteau et al., 2018; di Matteo et al., 2018). Future tests against catalogs with such a dataset could prove informative regarding correlations with extragalactic sources. Moreover, the ongoing upgrade of the Telescope Array, aimed at increasing the effective area of the observatory by a factor of 4 (Sagawa, 2016), will significantly reduce the contrast between the Northern and Southern exposures. Space-based observations with sufficient angular resolution could provide in the mid-term future a complementary approach to avoid UHECR blind spots over the celestial sphere.

Most current anisotropy studies exploit the arrival directions of UHECR events above a given (or scanned) energy threshold. This information could be supplemented by spectral and composition data to perform tomography of the UHECR production rate. Propagation of nuclei of different species affects the expected composition and spectrum as detected on Earth. Combined fits of the spectral and composition data show constraining power on the evolution of the density of sources at a fixed a luminosity (see e.g., Aab et al., 2017d). Constraints on composition are mostly inferred from fluorescence data, limited in statistics beyond few tens of EeV. With the upgrade of the Pierre Auger Observatory, the joint detection of showers with scintillators and water tanks will provide a composition-dependent observable with nearly 100% duty cycle (Aab et. al., 2016b). The selection of a “light” component, expected to be more localized than heavier nuclei, or the development and fit of models accounting for different propagation effects through diffuse photon and magnetic fields for different species could provide a clearer view on the population of sources (Alves Batista et al., 2015; Boulanger et al., 2018). Finally, interesting new approaches have emerged that aim to jointly model the UHECR spectrum and arrival directions, suggesting the possibility to associate a larger fraction of events to sources in catalog-based studies when accounting for the energy on an event-by-event basis (e.g., Farrar, 2008; Capel and Mortlock, 2019). These recent works suggest the possibility in the mid-term future to design analyses jointly accounting for the energy, composition, and arrival directions of UHECRs. This could for the first time enable a three-dimensional probe of the UHECR production rate, to be compared to the distribution of sources in the nearby Universe.

Assuming that sources of UHECRs also accelerate electrons radiating photons in a relativistic flow with speed β and bulk Lorentz factor Γ, the Hillas condition imposes a minimum photon luminosity L_γ which reads, under the assumption of equipartition between electrons and the magnetic field (Blandford, 2000; Lemoine and Waxman, 2009):

where the rigidity E/Z is currently estimated to be in the range 10¹⁸ − 10¹⁹ V beyond the ankle and where Γ²/β can range down to 10 for a mildly relativistic shock with β = 0.1 and Γ²/β = 100 either for Γ ~ 10, typical of blazar jets on pc scales, or for β ~ 10⁻², typical of starburst winds. UHECRs beyond the ankle could originate from sources up to about a Gpc. Then, the condition in Equation (4) corresponds to a minimum detectable flux for a full-sky electromagnetic survey at the level of $S_{\text{min}}=2\times 10^{-12}\text{ erg c}\text{m}{-}^2{\text{s}}^{-1}$, matching the current sensitivity limits of full-sky surveys from, e.g, Fermi-LAT in the gamma-ray band or WISE in the infrared¹. It thus appears that a census of potential UHECR sources beyond the ankle based on some electromagnetic full-sky surveys could be at hand. Two hurdles limit this statement. The first one is the sensitivity to extragalactic sources behind the Galactic plane, which acts as a strong foreground. The second one lies in the lack of full-sky spectroscopic surveys providing redshift information down to photometric sensitivity limits. While significant progress has been made in constraining the redshift distribution from electromagnetic surveys (see e.g., Bilicki et al., 2016; Cuoco et al., 2017 for recent contributions), further efforts may be needed to identify the best tracers of UHECR sources in a tomographic manner, accounting for incompleteness and possible contamination in every corner of the visible UHECR Universe.

3.2.4. Steady and Transient Sources

All known non-thermal sources are transient on some timescale. For UHECR sources, what defines whether a candidate object is classified as a transient or steady source is the ratio of the mean propagation timescale between sources to the source emission timescale, t_prop/t_emiss. For steady (transient) sources this ratio is < (>) 1. Both quantities, t_prop and t_emiss, are dependent on the UHECR energy.

The propagation timescale depends on the mode by which UHECRs propagate, which itself depends on the distance between sources and the UHECR scattering length. For a given source density, a CR energy can be found for which the distance between sources matches the cosmic-ray scattering length in the turbulent medium it is propagating through (Kotera and Lemoine, 2008b; Takami and Murase, 2012). Below this energy, cosmic rays propagate between sources on timescales significantly longer than than the ballistic propagation time.

The source emission timescale, t_emiss, is dictated by the collective timescale for particle acceleration, escape, and losses. On energetic grounds, only efficient Fermi acceleration (η ≲ 10) in sites associated to particular regions in AGN and GRB outflows satisfy the Hillas criterion in order to be considered as potential UHECR sources (see Equation 3). Much of what we know about these classes of astrophysical accelerators and their acceleration efficiency comes from the observation and analysis of their non-thermal emission (Murase and Takami, 2009).

For AGN, the longest timescale which may be associated to particle acceleration is the jet activity timescale, estimated to be of the order 300 Myr (Wykes et al., 2013). However, much shorter variability timescales are observed in the very-high energy gamma-ray (VHE,>100 GeV) emission of AGN. Studies of distant bright AGN sources over long epochs indicate that these objects release roughly an equal amount of power in logarithmic variability time bins over all epochs currently probed, from ~100 year down to daily timescales (Abdalla et al., 2017).

For GRBs, extensive efforts to detect VHE gamma rays have until recently failed to achieve a detection (Mirzoyan, 2019). Currently, the published record for the highest-energy emission observed is that seen at energies close to 100 GeV by Fermi-LAT, from the brightest GBM (in fluence) GRB event GRB130427A. The timescale for this emission was < 1,000 s. These results leave unclear whether GRBs operate as efficient particle accelerators (i.e., close to the Bohm limit), and on what timescale the acceleration takes place on. It therefore remains unclear whether these objects can be considered as viable UHECR sources.

Adopting a fiducial distance between sources of ~ 10 Mpc, a ballistic propagation time between sources of 30 Myr sets a lower limit to the actual propagation time. Adopting the 300 Myr AGN jet activity timescales as a fiducial value for t_emiss, only cosmic rays which diffusively scatter on a length scale >1 Mpc will contribute to the total flux as a steady-state contribution.

The above example demonstrates that the flux from any source class of a similar number density whose emission timescale is significantly shorter than a Myr will almost certainly be transient, and unable to achieve steady state. Furthermore, steady-state emission at low energies eventually becomes invariably unachievable for all source classes, once the diffusive sphere of cosmic rays around each source ceases to overlap with even neighboring sources—a phenomenon referred to as the magnetic horizon.

At high energies, energy losses during propagation affect whether the flux can achieve steady-state or not through a reduction of t_emiss. Should the reduction in this timescale lead to t_prop/t_emiss > 1, the flux at high energies will not be in the steady-state regime. Indeed, it is possible for a source class to only achieve steady-state emission in a finite energy range, both below and above which the cosmic-ray flux is only transient.

3.2.5. Origin of the End of the Cosmic-Ray Spectrum

The cut-off at the highest energies in the cosmic-ray spectrum has been established unambiguously recently, but the origin of this most prominent and significant feature is still a matter of debate. It has been tempting to identify the flux suppression with the long-predicted GZK-effect given its close coincidence to the expected threshold energy of about 6 · 10¹⁹ eV. Several fits of the end of the cosmic-ray spectrum with a propagated cosmic-ray composition consisting of a single element (p, He, N, Si, or Fe) at the source are shown in Figure 13. Different model lines are for different distances of the closest source and the source evolution was assumed to follow either the star formation rate from Robertson et al. (2015) or the AGN density from Ahlers et al. (2005). The cosmic-ray energy losses in the background photon fields were simulated with the CRPropa package (Alves Batista et al., 2016). In all panels, the spectral index at the source was a free parameter. Introducing additionally the maximum energy at the source as a free parameter (second and fourth panel) the measured spectra at Earth can be described well, no matter what is the composition at the source. However, “pure” GZK scenarios (first and third column) for which the maximum energy was fixed to 10²² eV and the flux suppression is only due to propagation effects do not fit the flux measured by the Pierre Auger Observatory well with the exception of Fe. For the TA measurement, which has larger statistical uncertainties and a higher energy scale, the GZK scenarios fit reasonably well.

A closer look at GZK scenarios is given in Figure 14, where a subset of fits from Figure 13 is shown. Here the sources are homogeneously distributed in the Universe following the star formation rate and emit pure beams of either protons, nitrogen, or iron nuclei with spectral indices γ chosen to fit the shape of the observed distribution. The source spectra are truncated exponentially at energies above 10²² eV so that the observed cut-off is caused by the GZK-effect. Also shown are the reconstructed compositions from Auger in terms of the mean mass 〈ln A〉 and its variance V(ln A). The lines show the results of the CRPropa simulations for the all-particle energy spectra (left), 〈ln A〉 (middle), and V(ln A) (right) for p, N, and Fe beams emitted from the sources. None of these simulations provides an acceptable description of the Auger and TA data. This is most obvious for the expected and observed compositions but also the simulated and observed all-particle spectra differ: in the Auger data sets, the suppression is below the GZK cut-off and in the TA datesets, all data points are above the GZK cut-off. The results question the interpretation of the flux suppression as caused solely by the GZK effect. While the energy spectra of Auger and TA could be made to agree with the GZK effect if the uncertainties of the energy scales are accounted for—14% in Auger and 21% in TA—the mass compositions of data and simulations for each experiment are totally different.

On the other hand, a simple astrophysical model of identical UHECR sources that accelerate nuclei through a rigidity-dependent mechanism provides a perfect description of the energy spectrum and mass composition above the ankle if the maximum rigidity is at about 10^18.8 V, the composition is dominated by intermediate mass nuclei, and the source spectra are harder (γ ≃ 1.6) than expected by the standard Fermi mechanism (Fang et al., 2013; Aab et al., 2017d; Wittkowski, 2018). In such a scenario, the flux suppression is a combination of propagation effects and the maximum energy at the source.

3.2.6. Magnetic Fields

Magnetic fields in scales comparable to and larger than the size of the Galaxy may affect the propagation of UHECRs. Little is known about extragalactic magnetic fields (EGMFs). The mechanisms whereby they originated are broadly divided into two classes, astrophysical and primordial. The latter postulates that fields in the present epoch result from the amplification of seed fields generated through a cosmological process in the early Universe, whereas in the former scenarios astrophysical processes such as feedback by active galaxies and stars would seed the intergalactic medium. Comprehensive reviews on cosmic magnetogenesis can be found in Durrer and Neronov (2013) and Subramanian (2016). Most of the Universe is filled with cosmic voids, whose magnetic fields are poorly constrained, ranging from 10⁻¹⁷G (Neronov and Vovk, 2010), estimated using gamma-ray-induced electromagnetic cascades, up to ~ 10⁻⁹G, from CMB measurements (Ade et al., 2016; Jedamzik and Saveliev, 2018). The lower bound, however, has been subject to much controversy (see e.g., Tiede et al., 2017; Broderick et al., 2018). In cosmic filaments, magnetic fields are B ~ nG, and in the center of galaxy clusters, B ~ μG; for reviews (see e.g., Vallee, 2011; Ryu et al., 2012).

The Galactic magnetic field (GMF) is understood better than EGMFs. Observationally driven models have been developed using polarized synchrotron maps, combined with Faraday rotation measurements (RMs) for the regular field and synchrotron intensity maps to derive the random field component. One of the most complete models was developed by Jansson and Farrar (2012a,b)—henceforth JF12. Yet, uncertainties are significant and improvements can be made. For instance, in Terral and Ferrière (2017) different models of the halo field were compared to the rotation measures. Different models for the disk and halo field were also studied in Unger and Farrar (2017) and in addition the uncertainties on the GMF due to the uncertainty of the synchrotron data and models for the thermal and cosmic-ray electrons were quantified. Theory-driven models of the GMF based on MHD simulations of structure formation can also provide complementary information and possibly improve the current picture (Pakmor et al., 2014; Pakmor et al., 2017; Marinacci et al., 2018) or alternatively dynamo-inspired models can be used to describe the large-scale GMF (Shukurov et al., 2018).

A detailed study of the deflection of UHECRs in the GMF model by JF12 has been performed (Farrar and Sutherland, 2017). By backtracking UHECRs at various energies to the edge of the Galaxy, the authors show that deflections for rigidities below 10 EV are large (≳ 90°). They also show that significant (de)magnification occurs for most of the rigidities studied (between 1 and 100 EV). The image patterns formed due to magnetic lensing have a considerable dependence on the ill-constrained turbulent component of the field (see also Battaner et al., 2011).

In Erdmann et al. (2016), the authors investigate UHECR deflections in the JF12 and Pshirkov et al. (2011) models. They point out that for cosmic-ray rigidities higher than 20 EV, these models lead to deflections compatible with each other, except near the Galactic disc (|b| ≤ 19.5°). A similar conclusion was reached by Unger and Farrar (2017) studying a larger ensemble of GMF models (19 variations of the JF12 model). If these different GMF models give a fair representation of the uncertainty of our knowledge of the deflections in the Galaxy, then it might be possible to correct the arrival directions of UHECRs for GMF-deflections above 20 EV.

Due to the lack of observationally derived models for the distribution of EGMFs in the local Universe, most studies of the kind have been done using cosmological simulations of structure formation. Early works by Sigl et al. (2003a,b) and Dolag et al. (2005) have reached conflicting conclusions regarding the role played by EGMFs on UHECR deflections. The latter concluded that deflections are small for E ≳ 40EeV, whereas the prospects for UHECR astronomy according to the former seem unfavorable. The origin of this discrepancy is related to the assumptions made, such as magnetogenesis mechanism (astrophysical or primordial), power spectrum of the seed magnetic field, local distribution of magnetic fields near the observer, among others. In Das et al. (2008) it has been argued that deflections would be <~ 5° in about a third of the sky. More recent works (Hackstein et al., 2016, 2018) considered both astrophysical and primordial magnetic field seeds. The authors attempt to cover the aforementioned uncertainties by studying the impact of different models of dynamo amplification and feedback by active galaxies, which may considerably change the distribution of magnetic fields. They confirm the predictions by Dolag et al., that UHECR deflections due to EGMFs are rather small. An extreme scenario with strong magnetic fields has been studied in Alves Batista et al. (2017), considering several magnetic power spectra for the seed fields. In this case, deflections of UHE protons with E ≳ 50 EeV are estimated to be <2° in about a quarter of the sky. For E ≳ 100EeV, nearly all protons would be deflected <~ 10°.

In addition to magnetic deflections, EGMFs induce energy and charge dependent time-delays as discussed in section 3.2.1. If the sources of UHECRs are transient, these time-delays are expected to produce an observable distortion of the arrival direction distribution with respect to that of steady sources (Kalli et al., 2011), a different energy dependence of the apparent UHECR source number density than steady sources (Takami and Murase, 2012), and spiky features at the highest energies of the UHECR spectrum from the brightest, most recent UHECR transients, that could help distinguish between steady and transient UHECR source populations (Miralda-Escude and Waxman, 1996).

Note, however, that magnetic power spectra which contain most of the energy at large scales would completely spoil UHECR astronomy in this case, though this is an unlikely scenario. Ultimately the deflection of UHECRs in EGMFs depends on the distribution of EGMFs in the local Universe. This can be understood in terms of volume filling factors, shown in Figure 15.

The effects of EGMFs on the spectrum and composition measured at Earth depend on the characteristic lengths involved. In the limit of a continuous source distribution, the propagation theorem states that the spectrum will have a universal form regardless of the modes of propagation (Aloisio and Berezinsky, 2004); this condition may not be easily realized, though. When the propagation time of a cosmic ray from its source to Earth is comparable to the age of the Universe, magnetic horizon effects suppress the lower energy region of the spectrum. In Mollerach and Roulet (2013) it was argued that this effect may play a role at E ≳ EeV for B ~ 1 nG. However, for realistic magnetic field distributions obtained from cosmological simulations, this effect may not be relevant at these energies depending on the source distribution and the distance to the nearest sources (Alves Batista and Sigl, 2014). Ultimately, it depends on the magnetic field distribution between Earth and the nearest sources.

A combined spectrum-composition fit of the Pierre Auger Observatory data including a particular model of EGMF has been presented in Wittkowski (2018). The results indicate a softening of the best-fit spectrum compared to the case where no magnetic field is assumed. This demonstrates the importance of understanding EGMFs in detail to improve phenomenological models.

Knowledge about the intervening magnetic fields is important to understand the origins of UHECRs. Conversely, UHECRs may also be used to constrain cosmic magnetic fields. A number of methods have been proposed for this purpose (Golup et al., 2009; Erdmann and Schiffer, 2010; Sutherland et al., 2012; Zimbres et al., 2014). In particular, Golup et al. (2009); Giacinti et al. (2010); Giacinti and Semikoz (2011); Zimbres et al. (2014); Oikonomou and Mostafa (2018) have proposed different searches for magnetically-aligned energy-ordered multiplets, which could be used to constrain the GMF, although none have been observed so far (Abreu et al., 2012).

Some attempts to constrain EGMFs using UHECRs have been made. Yuksel et al. (2012) have argued that at energies of ~100 EeV, a putative correlation of events observed by Auger with Cen A (Aab et al., 2015d) would imply an EGMF with strength B ≳ 20 nG, assuming a pure proton composition and the GMF model from Pshirkov et al. (2011). It is also possible to constrain properties of EGMFs other than the magnetic field strength. For instance, in Alves Batista and Saveliev (2019) a method has been proposed to infer the helicity—a topological quantity related to the degree of twisting and linkage of magnetic field lines — of cosmic magnetic fields. Although helicity is often neglected in UHECR studies, it has been proven to leave imprints in the large-scale distribution of UHECRs.

Given the considerable uncertainties in the GMF and the inconclusive results on the effects of EGMFs on cosmic-ray propagation due to the model dependencies, the unambiguous identification of individual UHECRs will require better constraints on EGMFs and improved models of the GMF. Therefore, until we have a more accurate description of GMFs and unless EGMFs are small or better understood, charged-particle astronomy will remain challenging if cosmic rays at the highest energies are heavy nuclei until new data on the RMs of Galactic pulsars and Faraday tomography [e.g., from LOFAR and SKA (Beck, 2008)]. Large-scale surveys of starlight polarization (Magkos and Pavlidou, 2019) will soon allow for a better three-dimensional reconstruction of the GMF.

3.3. Particle Physics

3.3.1. Hadronic Interactions at Ultrahigh Energies

Even though the current generation of hadronic interaction models gives a good description of many properties of air showers, we are far from having reached a satisfactory level in the quality and reliability of modeling extensive air showers.

First of all, there is the muon discrepancy that is still not understood. An excess in the muon number of ~ 30% relative to simulation predictions is not accounted for (Aab et al., 2015b; Aab et al., 2016a). This is the most important question to be addressed in model development and forthcoming fixed-target and collider measurements. The difficulty is that so far no “smoking gun” signature has been found that might indicate in what direction to search. It will be the task of further air-shower measurements to characterize the muon discrepancy in terms of the energy spectrum, production depth, and lateral distribution. Also going beyond measuring mean values will be important, as muon fluctuations give a handle on features of the first few interactions in a shower (Cazon et al., 2018a,b). Conventional explanations of the muon discrepancy have to be developed to find out whether we indeed have to assume that the muon excess is related to physics beyond the Standard Model of particle physics.

Secondly, the accuracy of the predictions of hadronic interaction models has to be improved to reduce the systematic uncertainties of composition measurements. For example, the tension between the mean depth of shower maximum and the shower-by-shower fluctuations of X_max imply an almost mono-elemental mass composition. The astrophysical implications of such an unexpected change of mass composition with energy, without reaching a truly mixed composition, are severe (see e.g., Aab et al., 2017d). Therefore, it is of prime interest to figure out whether this apparent tension is not an artifact of inaccurate expectations for the mean X_max. Moreover, the comparison of the electromagnetic and muonic shower components at surface level offers an even larger composition sensitivity than that provided by the depth of shower maximum. Currently, this sensitivity cannot be fully exploited. Muon-based mass measurements systematically lead to heavier mass compositions than X_max-based analyses (Kampert and Unger, 2012).

Over the last decades, a rich dataset on proton-proton and proton-antiproton interactions at high energy has been accumulated in various collider experiments. In an air shower, most of the interactions are initiated by pions and kaons, except the first one. There is a severe lack of pion-proton and kaon-proton data that are needed for improving our understanding of hadronic interactions and for tuning interaction models. Of prime importance are the measurement of the pion-proton, kaon-proton, pion-light-nuclei, and kaon-light-nuclei cross sections and the corresponding distributions of leading secondary particles. Taking data at LHC and selecting events in which a beam proton becomes a neutron by emitting a π⁺ is one possibility to study pion interactions at energies not accessible in fixed target experiments. Similar measurements have been done at HERA (Khoze et al., 2006). Such measurements could help determine the secondary particle distributions, but the interaction cross section of pions and kaons can only be measured in fixed-target experiments. Additional information could be obtained from air showers by studying the muon production depth (Ostapchenko and Bleicher, 2016).

Similarly, there is almost no data available on particle production with light nuclei in the mass range of air (〈A〉 ≈ 14.45). The current understanding of nuclear processes is too limited to be able to reliably predict the secondary particle distributions in proton-air interactions if they were known for proton-proton and proton-neutron interactions. Also, the measurements of heavy-ion interactions, such as p-Pb and Pb-Pb, cannot be transferred to light nuclei with the needed accuracy. Taking data of proton-oxygen at LHC is technically possible and would be a key measurement for improving air shower predictions.

And, last but not least, particle detectors covering the forward direction would help significantly to reduce the needed extrapolation of collider measurements to phase space regions of relevance to air showers. Understanding the scaling of the forward particle distributions at collider energies is the key to extrapolating to higher interaction energies. LHCf (Adriani et al., 2018) is a good example that such detectors can be built even though there are large technical challenges and limitations.

Although it can be expected that progress in understanding hadronic multiparticle production will be mainly driven by experimental results for the next years, efforts to develop a more consistent theoretical framework will be equally important. The transition between soft and hard processes (i.e., processes with small and large momentum transfer) is not understood at all. Applying Regge parameterizations for soft processes is common practice, while hard processes are treated within the QCD-improved parton model. Closely related to this transition between two regimes is the question of non-linear effects, or even parton density saturation, expected at very high parton densities. LHC data show that proton-proton interactions of high multiplicity exhibit features previously only seen in heavy ion collisions (see e.g., Khachatryan et al., 2010; Aad et al., 2016c; Adam et al., 2017). Effects related to high parton densities will have to also be considered in hadronic interaction models for air showers.

3.3.2. Physics Beyond the Standard Model

The center-of-mass energy of an UHECR of energy E, in the lab frame, interacting with a low-energy particle of energy ϵ is given by

Thus, UHECR propagation in the cosmic radiation backgrounds, which have energy ϵ≪ GeV, cannot probe Lorentz-invariant physics beyond the Standard Model, but it can still probe invariance violations under Lorentz boosts. In addition, the development of air showers can be influenced both by Lorentz-invariant new physics acting on primaries or secondaries with energy E ≳ 2 × 10¹⁷eV and targets of mass ϵ ≳ m_N, and also by Lorentz-invariance violations involving large Lorentz boosts.

Air-shower physics. Observed air showers show an excess of muons compared to predictions of standard hadronic interaction models above ≃ 10¹⁶ eV in the KASCADE-Grande experiment, which indicates a longer-than-expected muon attenuation length (Apel et al., 2017). The Pierre Auger Observatory also sees a muon excess by a factor ≃1.5 (Aab et al., 2016a). Since a large number of muons is observed, this cannot be a statistical effect but rather points to shortcomings in the models. If this muon excess cannot be explained by improved hadronic event generators within the Standard Model, new physics could qualitatively play a role in the following way.

The muon abundance is roughly proportional to the energy fraction going into hadrons. Since over one hadronic interaction depth of about $X_h(E)\simeq \left[88-9log(E/\text{EeV})\right]\text{g}{\text{ cm}}^{-2}$ a fraction $f_{{\pi }^0}$ — the branching ratio of hadronic interactions into neutral pions—of the hadronic energy is converted into electromagnetic energy, after n hadronic generations the hadronic energy fraction of the shower is proportional to ${(1-f_{{\pi }^0})}^n$. Therefore, a larger muon number could be caused by decreasing the number of generations n, decreasing $f_{{\pi }^0}$, or both. However, the number of generations is well constrained by detailed measurements of X_max. Therefore, assuming neutral pions still decay quasi-instantaneously, the most likely explanation of the observed large muon number is a significant decrease in the fraction of energy going into neutral pions, $f_{{\pi }^0}$. For example, it has been suggested (Farrar and Allen, 2013a,b) that if chiral symmetry is restored above a certain center-of-mass energy, pions may become much heavier and their production may be suppressed in favor of baryon-anti-baryon production. This would put more energy into the hadronic channel compared to the electromagnetic channel, thus producing more muons. A similar effect could be achieved by the production of a fireball consisting of deconfined quarks and gluons (Anchordoqui et al., 2017).

Alternatively, if high-energy neutral pions were stable or had a decay rate smaller than their interaction rate in the atmosphere, for example, due to Lorentz symmetry violation at very high Lorentz factors, then their energy could contribute to increasing the energy fraction going into the hadronic channel and thus the muon signal. Generally speaking, an increase of the fraction of air-shower energy in the hadronic channel would likely be a hint for new physics. One could imagine such effects to have energy thresholds, so one could also search for comparatively large increases of muon number over a small primary energy range.

Lorentz-invariance violation. Lorentz invariance violation (LIV) can be induced by non-renormalizable operators that conserve gauge invariance but break parts of the Poincaré group (Colladay and Kostelecky, 1998). For example, it has been shown that in quantum electrodynamics the most general non-renormalizable dimension-five CPT−odd operator that is quadratic in the fields and preserves rotation and gauge invariance, but is not invariant under Lorentz boosts, can be written as (Myers and Pospelov, 2003)

Here, ξ, χ₁ and χ₂ are dimensionless constants, M_Pl is the Planck mass, u^μ is a constant time-like four-vector which corresponds to a preferred Lorentz frame such as the cosmic microwave background rest frame and ${\tilde{F}}^{\mu \nu }$ is the dual electromagnetic field strength tensor.

Operators such as Equation (6) can manifest through modifications of dispersion relations for particles of energy E, momentum p, and mass m, by terms that are suppressed by a power n of the Planck mass M_Pl(Amelino-Camelia et al., 2005; Christian, 2005; Diaz, 2014). The dispersion relation for left- and right-handed photons or fermions can be written as

Here, n = d − 4 for a d-dimensional operator and the dimensionless numbers η_± refer to positive and negative helicity states, respectively. In general, in effective field theory one has ${\eta }_{+}={(-1)}^n{\eta }_{-}$. For example, Equation (6) implies n = 1 and η_± = ±ξ for right- or left-circularly polarized photons, respectively, and η_± = 2(χ₁ ± χ₂) for positive and negative electron helicity, respectively. For renormalizable LIV terms, d ≤ 4 and n is negative in Equation (7).

Dispersion relations of the form of Equation (7) can modify both the free propagation of particles and the kinematics and thresholds of interactions (Amelino-Camelia and Piran, 2001; Dubovsky and Tinyakov, 2002; Moffat, 2003; Gagnon and Moore, 2004; Stecker and Scully, 2005; Galaverni and Sigl, 2008; Bi et al., 2009; Maccione et al., 2009; Scully and Stecker, 2009; Cowsik et al., 2012). Kinematics are typically modified when the LIV terms become comparable to the particle rest mass, i.e., when the particle energy is larger than a critical energy E_cr,

Therefore, the larger the particle mass, the higher the energy at which LIV effects become relevant.

In the relativistic limit, to first order in m² and η_±, the group velocity corresponding to Equation (7) is

For positive η this would lead to superluminal motion for E > E_cr. Particles with $v_{\text{gr}}^{\pm }$ tend to emit vacuum Cherenkov radiation, similar to the motion of an ultra-relativistic charge in a medium with index of refraction larger than one, and would lose energy rapidly (Cohen and Glashow, 2011). Therefore, observing a particle of energy E implies E_cr ≳ E which, from Equation (8), leads to an upper bound on η. For n = 1, the observation of EeV protons places a stringent limit on LIV, of

If LIV exists and its effect is non-negligible, the corresponding parameter η should naturally be of order 1; constraining it to values much smaller than unity (Gagnon and Moore, 2004; Galaverni and Sigl, 2008) would suggest that the corresponding LIV does not exist; however, see Aloisio et al. (2014). Additionally, Equation (9) leads to energy-dependent delays in the propagation time from the sources to Earth. Thus, strong constraints on LIV may be placed by the detection of UHE photons from local sources (Murase, 2009).

New physics in UHE neutrinos. High-energy astrophysical neutrinos, with TeV–PeV energies, recently discovered, have opened up a new regime to test for new physics (Anchordoqui et al., 2005; Gonzalez-Garcia et al., 2005; Ahlers et al., 2018). They are unparalleled in two key features: they have the highest neutrino energies detected—so they can probe effects at new energy scales— and they travel over the longest baselines—so tiny new-physics effects could accumulate en route to us, and reach detectable levels. Cosmogenic neutrinos, with EeV energies, when discovered, could extend the reach of these tests.

IceCube astrophysical neutrinos have been used to measure the neutrino-nucleon cross section in the TeV–PeV range for the first time (Aartsen et al., 2017; Bustamante and Connolly, 2019). It was found to be compatible with high-precision Standard Model predictions based on collider data (Cooper-Sarkar et al., 2011), though there is still room for small deviations due to new physics. Cosmogenic neutrinos could be used to measure the cross section at the EeV scale for the first time, test strong dynamics more deeply than colliders (Anchordoqui et al., 2019; Bertone et al., 2019), and search for new physics at an even higher energy scale.

Numerous new-physics models have effects that are proportional to some power of the neutrino energy E and to the propagated distance L, i.e., they grow as $~{\kappa }_n{E}^{n}L$, where the energy dependence n and the proportionality constant κ_n are model-dependent. For instance, for neutrino decay (Chikashige et al., 1980; Gelmini et al., 1982; Tomas et al., 2001), n = −1; for CPT-odd Lorentz violation (Colladay and Kostelecky, 1997; Coleman and Glashow, 1999; Barger et al., 2000) or coupling to a torsion field (De Sabbata and Gasperini, 1981), n = 0; and for CPT-even Lorentz violation (Coleman and Glashow, 1997; Glashow et al., 1997) or violation of the equivalence principle (Gasperini, 1988, 1989; Halprin et al., 1996; Adunas et al., 2001), n = 1. An experiment that sees neutrinos of energy E coming from sources located at a distance L is, in principle, able to probe new physics with sensitivities of ${\kappa }_n~4·10^{-50}{(E/\mathrm{\text{EeV}})}^{-n}{(L/\mathrm{\text{Gpc}})}^{-1}{\mathrm{\text{EeV}}}^{1-n}$, a significant improvement over current limits of ${\kappa }_0\lesssim 10^{-29}$ PeV and ${\kappa }_1\lesssim 10^{-33}$ (Abbasi et al., 2010; Abe et al., 2015).

New physics of different types can affect all neutrino observables: the energy spectrum (see e.g., Baerwald et al., 2012; Blum et al., 2014; Ibe and Kaneta, 2014; Ioka and Murase, 2014; Ng and Beacom, 2014; Kopp et al., 2015; Bustamante et al., 2017), distribution of arrival directions (see e.g., Davis and Silk, 2015; Argüelles et al., 2017), and the flavor composition, i.e., the proportion of each neutrino flavor to the total in the incoming flux (see e.g., Argüelles et al., 2015; Bustamante et al., 2015; Shoemaker and Murase, 2016; Rasmussen et al., 2017).

At high and ultra-high energies, there are a few challenges to detecting new physics in high-energy cosmic neutrinos:

• New-physics effects might be sub-dominant; in this case, their discovery is contingent on detecting a large enough number of events and on accurately reconstructing key properties of events, like energy and arrival direction;

• When extracting fundamental neutrino properties from the data, one must factor in astrophysical uncertainties (e.g., shape of the energy spectrum, redshift evolution of the number density of sources, etc.), which can be significant;

• Flavor is a difficult property to measure in neutrino telescopes (Aartsen et al., 2015a,b; Vincent et al., 2016; Bustamante and Ahlers, 2019); improved methods of flavor identification might be needed to fully exploit flavor in studying fundamental neutrino physics (e.g., Wang et al., 2013; Li et al., 2016).

These challenges are likely surmountable. High- and ultra-high-energy observatories are in a unique position to perform powerful tests of neutrino physics, complementing and expanding tests performed by experiments with lower energies and shorter propagation baselines.

4. Discussion

4.1. The Current Status and Perspectives of Earth-Based UHECR Detectors

The most promising step in the activity of ground UHECR detection is the upcoming upgrade of the Pierre Auger Observatory (Aab et. al., 2016b). It consists in the improvement of the surface detector (SD), namely water-Cherenkov tanks, by equipping each tank with solid-state scintillator plates on top. This configuration allows one to improve the sensitivity to the CR mass composition by simultaneous measurements of electrons and muons passing through both detectors. Although the total aperture of the surface detector will not increase, the fraction of events with a reliable reconstruction of the mass composition will be larger; day-time SD data as well as night-time data recorded by the fluorescence detectors (FD) will be cross-checked with new scintillators, thus improving the quality of hybrid SD+FD events.

Another upgrade of the Pierre Auger Observatory, recently confirmed, is the equipment of the SD with radio antennas. Contrary to the existing AERA detector (Aab et al., 2018b), which densely covers only a small part of the observatory, the new detector will feature a sparser layout and cover the full area of the observatory. In recent years, different antenna types for air-shower detection (Abreu et al., 2012) were investigated by a number of experiments. Based on these studies, the loop antenna, which was successfully exploited at the Tunka-Rex experiment (Bezyazeekov et al., 2015, 2018), was selected for the Pierre Auger upgrade. Joint operation of particle and radio detectors decreases the systematic uncertainty of energy and mass composition reconstruction, since radio detection allows one to reconstruct the calorimetric energy of the electromagnetic part of the air shower as well as the depth of shower maximum.

Also TA has recently started to be upgraded (Kido, 2018). Once the upgrade is complete, the array, TA × 4, will consist of three times more surface detectors than TA, similar to the original ones (two solid-state scintillators separated by a metal plate). The upgraded detector will cover an area of about 3 000 km², with the new scintillators two times sparser than the old ones. Additional FD will be built for hybrid operation with the extended array. The aperture of TA × 4 will facilitate the study of anisotropies at ultrahigh energy in the Northern Hemisphere and aid the comparison of the spectrum in the two hemispheres at the highest energies.

The detection of ultrahigh energy cosmic rays is also included in the scientific program of GRAND, the most ambitious ground-based experiment proposed so far (Alvarez-Muniz et al., 2018). Since the detector will consist of antennas tuned for the detection of very inclined events, its exposure overlaps both with TA and Auger. Due to the unprecedented exposure of GRAND in its envisaged final configuration (200 000 km²), it will be possible to detect about 32 000 cosmic-ray events with E > 10^19.5 eV in five years. Since GRAND exploits the radio technique for air-shower detection, cosmic-ray properties will be studied via measurements of the calorimetric energy of the air-shower and should achieve good X_max resolution.

4.2. The Current Status and Perspectives of Space Experiments to Study UHECRs

J. Linsley and R. Benson were the first to propose measurements of the fluorescent radiation of EAS using a UV telescope on-board a satellite (Benson and Linsley, 1981). Y. Takahashi, later proposed the idea of using wide-angle optics and CCD readout in the MASS concept (Takahashi, 1995). A space-based detector for UHECR research has the advantage of a much larger exposure and uniform coverage of the celestial sphere. This idea has been developed in a number of projects. In the late 1990s, the Airwatch concept was developed by J. Linsley, B. Scarsi, Y. Takahashi and others based on Fresnel optics. They later collaborated with a team from Utah/GSFC who separately developed the OWL/Crystal Eye idea to propose OWL-Airwatch (Streitmatter, 1998), a concept for a 2-spacecraft mission. The OWL concept later moved to Schmidt telescopes and into the final OWL study.

The original Airwatch concept, developed into the Extreme Universe Space Observatory (EUSO) (Catalano et al., 2001). This was the start of the JEM-EUSO program which originally took its name from the Japanese Experiment Module (JEM) but currently stands for Joint Experiment Missions. In the JEM-EUSO Collaboration, a large Fresnel lens telescope was developed (Adams et al., 2015). In Russia, detectors that use concentrator mirrors for collecting fluorescence light, TUS (Klimov et al., 2017) and KLYPVE (Panasyuk et al., 2016), were proposed and developed.

The TUS experiment was the first orbital detector of UHECRs. It was launched on board the Moscow State University (MSU) satellite “Lomonosov” (Klimov et al., 2017) on 28 April 2016. TUS is a UV telescope looking downward into the atmosphere in the nadir direction. It consists of two main parts: a modular Fresnel mirror-concentrator and 256 photomultiplier tubes (PMTs) arranged in a 16 × 16 photodetector located in the focal plane of the mirror. The overall field of view (FOV) of the detector is 4.5° × 4.5°. During 1.5 years of operation in EAS mode about 200 000 events of various types were measured during the night-time part of the orbit. The events differ in the spatial dynamics and temporal structure of their waveforms. Some EAS candidates have been registered.

Another, much larger space instrument, KLYPVE, is being developed in close cooperation with the JEM-EUSO Collaboration and is known as KLYPVE-EUSO (K-EUSO) (Panasyuk et al., 2016). To fulfill the requirements of the K-EUSO experiment, a Schmidt UV telescope covering a FOV of 40° with an entrance pupil diameter of 2.5 m, and a 4 m diameter mirror was developed. The baseline variant consists of a spherical mirror, a corrector plate and a spherical focal surface concentric with the mirror, placing the aperture stop on the frontal surface of the corrector plate. Even though the expected statistics of UHECR events will not exceed those of upcoming on ground installations (see Figure 16), with the current design, the K-EUSO instrument can perform the first all-sky observation of UHECRs, in order to establish whether the particle fluxes of the two hemispheres are different.

K-EUSO will measure about 140 UHECR events in the Northern hemisphere and 30 events in the Southern hemisphere at E > 57 EeV in one year of observations if the difference of the TA and Auger spectra is due to different fluxes. In contrast, the numbers of events from both hemispheres are expected to be nearly equal if the flux is isotropic. Also, K-EUSO data will allow for a full-sky search for UHECR anisotropy to independently confirm or rule out the presence of hotspots in the Northern and Southern hemispheres.

The project concept of OWL, based on the simultaneous detection of UHECRs by UV telescopes placed on two satellites, was recently developed in the POEMMA project (Olinto et al., 2018). This project, based on the use of Schmidt optics with 45° FOV and a large photodetector camera, can become a space instrument of record characteristics and surpass in terms of exposure the ground-based Auger and TA installations (see Figure 16).

4.3. The Current Status and Perspectives of UHE Neutrino Experiments

Currently the UHE neutrino flux is best confined by the IceCube Observatory (Aartsen et al., 2018) and the Auger Observatory (Bellido, 2018) at the level of ~ 3 × 10⁻⁸ GeVcm⁻²s⁻¹sr⁻¹ around EeV (all-flavor). Figure 17 summarizes the sensitivity of current and proposed experiments that target EeV neutrinos. The Askaryan Radio Array (ARA) (Allison et al., 2012, 2016) and ARIANNA (Barwick et al., 2017; Nelles, 2018) are in-ice radio arrays which detect UHE neutrinos via the Askaryan effect. As an alternative to the expensive ice-Cherenkov technique the three experiments equipped with radio antennas are located in Antarctica and optimized for UHE neutrino detection, namely two in-ice arrays, the Askaryan Radio Array (ARA) (Allison et al., 2012, 2016) and ARIANNA (Barwick et al., 2017; Nelles, 2018), and a balloon-borne interferometer ANITA (Gorham et al., 2018a; Gorham et al., 2018b). The propose GRAND (Alvarez-Muniz et al., 2018) will use large arrays of cost-effective radio antennas to detect particle cascades produced in media and air by UHE tau neutrinos. POEMMA (Olinto et al., 2018) will also detect tau neutrinos, by observing the Cherenkov radiation produced by upward-going tau decays (Neronov et al., 2017). Trinity (Otte, 2018), an Earth-based imaging telescope experiment, will detect air showers induced by taus or tau neutrinos by observing the Cherenkov or fluorescence light produced by the EAS.

5. Outlook

Despite revolutionary progress, some critical, long-standing questions in the field of UHECRs remain unanswered, or only answered partially: What are the sources of UHECRs? What is the mass composition of UHECRs at the highest energies? What mechanism accelerates CRs beyond PeV energies? What is the flux of secondary messengers—neutrinos, gamma rays—associated with UHECRs, and what can we infer from them about UHECR sources?

Observations performed by current and planned ultrahigh-energy facilities have an opportunity to give definite answers to these questions. Yet, to fulfill this potential, it is necessary to undertake a number of essential steps toward experimental and theoretical progress. Below, we list what we believe are the most important of these. This list is, of course, non-exhaustive and only expresses our views.

• UHECR composition: Precise measurement of the UHECR mass composition near the end of the spectrum is hindered by uncertainties in models of hadronic interaction, uncertainties in measuring X_max, and small statistics. The latter issue will be addressed by upgraded configurations of current facilities and larger, next-generation facilities.

Action item 1: Craft a program of accelerator measurements of cross sections and multiplicities to reduce uncertainties in models of hadronic interaction.

Action item 2: Explore experimental methods to infer the composition with precision comparable to that of fluorescence detectors and at a duty cycle of 100%.

• Identification of UHECR sources: Sources of UHECRs can be searched for either by self-correlation of arrival directions or by searching for positional correlations with sources catalogs.

Action item 3: Thorough studies of the effects Galactic and extragalactic magnetic fields.

Action item 4: When looking for correlations with source catalogs, generate catalogs providing a tomographic mapping of possible UHECR sources, for viable source populations, accounting for incompleteness and bias, up to the GZK radius at the energy of the ankle, possibly down to the minimum luminosity imposed by theoretical criteria.

• Particle acceleration: An accurate understanding of particle acceleration in astrophysical sources could help to interpret the transition from Galactic to extragalactic origin of cosmic rays and the shape of the UHECR spectrum at the highest energies, and would influence predictions of spectra of cosmogenic secondaries.

Action item 5: Perform detailed studies of particle acceleration in collisionless shocks and magnetic reconnection under conditions as close to those in real sources, either via simulations or in the lab, when possible.

• Muon excess problem: The mismatch between the number of muons with energies above 10^9.5 GeV predicted by shower models and the number detected points to further problems with the hadronic interaction models. Independent measurements of the electromagnetic and muon component of air showers could help solve this issue.

Action item 6: Favor the construction of air-shower facilities with separate electromagnetic and muon detectors.

• Updated prediction of cosmogenic neutrinos: Recent recalculations of the predicted flux of cosmogenic neutrinos, fitting the latest UHECR data, have resulted in fluxes significantly lower than before. Yet, the uncertainties in the prediction are large. This represents a problem in planning for the next-generation of UHE neutrino detectors.

Action item 7: Generate predictions of the cosmogenic neutrino flux by scanning across all of the available parameter space of UHECR model parameters—including uncertainties in magnetic fields and hadronic interaction models—in order to fully characterize the uncertainties.

• Updated predictions of UHE neutrinos from point sources: Because the flux of cosmogenic neutrinos might be tiny, UHE neutrinos from point sources might be detected first in next-generation neutrino telescopes. However, the literature on models of UHE neutrinos is outdated or lacking.

Action item 8: Generate updated predictions of the emission of UHE neutrinos from point sources, steady-state and transient.

• Global studies: A complete picture of the high-energy Universe needs to account for all messengers

Action item 9: Assess the validity of UHECR models by considering the full UHECR, neutrino, and photon data, from as many experiments as possible; avoid picking and choosing observables and experiments.

• Open data policies: For progress to be faster, the community should have access to detected events in UHE facilities, in a usable, non-raw form.

Action item 10: Existing and future facilities should have an open data policy, including software analysis tools when possible.

More than five decades of experimental and theoretical progress in the field of UHECRs will soon be compounded on by upgrades of Auger and TA, and by a suite of potential next-generation detectors. On one hand, thanks to these, in the next 5–10 years the increased statistics of UHECRs alone will refine the measurement of the energy spectrum, mass composition, and anisotropies to the point where several of the open questions above could already be answered. Additional improvements in analysis techniques will only enhance these prospects. On the other hand, upcoming detectors will potentially trigger a transformative change in the field: for the first time, we could reach the sensitivity needed to discover even tiny fluxes of cosmogenic neutrinos and gamma rays. Opening up the full breadth of UHE multi-messenger observables could answer most of the remaining open questions, and finally, provide a complete picture of the Universe at the highest energies.

Author Contributions

RB, JB, MB, RE, KF, K-HK, DK, GS, FO, MP, AT, and MU contributed to the original material and writing of the manuscript. FO and KF coordinated this review. All authors contributed to the discussions at MIAPP, read the manuscript and provided critical feedback.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge the support of the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence Origin and Structure of the Universe, where this work was initiated. We thank the organizers of the The High-Energy Universe: Gamma-Ray, Neutrino, and cosmic-ray astronomy MIAPP Program, Francis Halzen, Angela Olinto, Elisa Resconi, and Paolo Padovani, for the very fruitful workshop.

RB is supported by grant #2017/12828-4, São Paulo Research Foundation (FAPESP). MB is supported by the Danmarks Grundforskningsfond Grant 1041811001 and Villum Fonden project no. 13164. RE, K-HK, GS, and MU are supported by the Bundesministerium für Bildung und Forschung (BMBF) and the Deutsche Forschungsgemeinschaft (DFG). KF acknowledges support from the Einstein Fellowship from the NASA Hubble Fellowship Program. The work of KM is supported by Alfred P. Sloan Foundation and NSF grant No. PHY-1620777. FO is supported by the Deutsche Forschungsgemeinschaft through grant SFB 1258 Neutrinos and Dark Matter in Astro- and Particle Physics.

Footnotes

References