Probing UHECR sources - Constraints from cosmic-ray measurements

TL;DR

This paper reviews how recent cosmic-ray measurements and modeling advancements constrain the origins of ultra-high-energy cosmic rays, highlighting the role of extragalactic sources, composition, and anisotropy in understanding their sources.

Contribution

It provides a comprehensive review of current constraints on UHECR sources based on recent observational data and theoretical modeling advancements.

Findings

Extragalactic sources emitting intermediate-mass nuclei explain spectrum and composition.

Anisotropy measurements constrain source density and distribution.

Open questions remain about source similarity and highest energy event origins.

Abstract

Ultra-high-energy cosmic rays (UHECRs) are the most energetic particles known - and yet their origin is still an open question. However, with the precision and accumulated statistics of the Pierre Auger Observatory and the Telescope Array, in combination with advancements in theory and modeling - e.g. of the Galactic magnetic field - it is now possible to set solid constraints on the sources of UHECRs. The spectrum and composition measurements above the ankle can be well described by a population of extragalactic, homogeneously distributed sources emitting mostly intermediate-mass nuclei. Using additionally the observed anisotropy in the arrival directions, namely the large-scale dipole >8 EeV as well as smaller-scale warmspots at higher energies, even more powerful constraints on the density and distribution of sources can be placed. Yet, open questions remain - like the striking similarity of the sources that is necessary to describe the rather pure mass composition above the ankle, or the origin of the highest energy events whose tracked back directions point towards voids. The current findings and possible interpretation of UHECR data will be presented in this review.

Equations

$$Q_{A}(E)=Q_{0,A}\ \Big{(}\frac{E}{E_{0}}\Big{)}^{-\gamma}\ f_{\mathrm{cut}}(E,...).$$

Full text

1 Introduction

Although ultra-high-energy cosmic rays (UHECRs), nuclei with energies $>10^{18}\,\mathrm{eV}=1\,\mathrm{EeV}$, were first detected decades ago, their origin still remains an open question. What kind of powerful objects can accelerate particles to such energies, and how exactly the acceleration works, is still unclear. Through theoretical considerations, several constraints can be placed on UHECR sources, for example the Hillas criterion [1], which gives a necessary geometrical condition on the size and magnetic field of the accelerator for the case of electromagnetic acceleration. Additionally, the sources must have sufficient power to generate particles of the observed energies [2, 3, 4]. Another criterion can be set by requiring the sources to supply the observed UHECR intensity, implying that the product of their number density and luminosity must be large enough. For more information see e.g. these reviews [5, 6, 7, 8]. Different types of possible steady or transient UHECR sources have been proposed and discussed in the literature, for example active galactic nuclei (AGNs), gamma-ray bursts (GRBs), starburst galaxies (SBGs), galaxy clusters, tidal disruption events (TDEs), binary neutron star (BNS) mergers, and more, see e.g. [6, 9]. Several of these candidates fulfill the theoretical conditions, so definite conclusions regarding the sources of UHECRs must rely on measured data.

In the last years, a large increase in statistics has been achieved thanks to the two current biggest observatories for UHECRs measurements, the Pierre Auger Observatory [10] (Auger in the following) in the southern hemisphere, and the Telescope Array [11] (TA) in the northern one. Due to their massive areas and extensive measurement times since 2004 and 2008, respectively, precise measurements of the energy spectrum and arrival directions of UHECRs are now available. Additionally, the mass composition can be constrained by means of the depth of shower maximum or the amount of muons in the shower [12]. Still, the interpretation of the data remains a delicate task. This is in part due to the influence of magnetic fields that deflect UHECRs on their way from the sources to Earth. Recently, significant progress has been made regarding models of the Galactic magnetic field (GMF) [13, 14], enabling far more precise predictions of the directional deflections of UHECRs. Nevertheless, the influence of the extragalactic magnetic field (EGMF) remains quite unknown (for a review see [15]). An upper bound on the field strength in voids from observations of the cosmic microwave background is $\mathcal{O}(1\,\mathrm{nG})$ [16], but it is only applicable for a primordial origin of the EGMF. A comparable upper limit, independent of the EGMF origin, can however be placed using rotational measures [17]. Next to the influence of the GMF and EGMF, also effects such as the interactions with background photon fields have to be considered when interpreting UHECR data. In recent years, models for the origin of UHECRs have become more and more sophisticated, allowing for better conclusions on UHECR sources, but also opening up even more questions.

In this review, recent progress regarding the interpretation of UHECR data will be discussed. For that, first the UHECR energy spectrum and composition data as well as their possible interpretation using models are presented in sec. 2, and important common conclusions about UHECR sources will be identified. In sec. 3, further constraints from arrival direction data will be discussed, starting first with a summary of the effects of the Galactic magnetic field in sec. 3.1. Possible interpretations of the arrival direction data will be explored, moving from large-scale anisotropies in sec. 3.2 over intermediate scales in sec. 3.3 to the highest energy events in sec. 3.4. In the end, a conclusion and outlook to the future will be given in sec. 4.

2 Characteristics of UHECR sources - constraints from the measured spectrum and composition

One of the key measurement of UHECRs is the energy spectrum, where several features have been significantly identified at different energies: the ankle at around $10^{18.7}$ eV, the instep or shoulder at around $10^{19.1}$ eV, and a suppression of the flux beyond $10^{19.7}$ eV. While the spectra measured in the southern hemisphere by Auger [18] and the northern hemisphere by TA [19] agree about the general features and align well at lower energies when considering systematic uncertainties, there is some discrepancy at higher energies [20] whose origin is still under debate.

The measurement of the depth of the shower maximum $X_{\mathrm{max}}$ allows for conclusions on the primary mass of UHECRs. The Auger fluorescence detector data shows a transition from a mixture of lighter elements below the ankle towards continuously heavier mass at higher energies [21, 12]. The second moment of $X_{\mathrm{max}}$ decreases continuously $\gtrsim 5\,\mathrm{EeV}$, indicating that the amount of mixing between elements decreases with energy and the composition becomes progressively purer. These trends were recently confirmed by surface detector data using a neural-network based analysis [22] which allowed for a 10-fold increase of statistics and an extension of the $X_{\mathrm{max}}$ moments to higher energies. Also, it was shown that the Auger data is better described if the predicted $X_{\mathrm{max}}$ by hadronic interaction models is shifted toward deeper values, implying an even heavier mass composition than previously inferred [23]. While the TA data is often interpreted as light even above the ankle [24], it was recently shown that it is compatible within uncertainties with both an all-proton as well as the heavier Auger mix interpretation [25]. The arrival directions measured at TA show a high level of isotropy implying a heavy mass composition, at least at the highest energies $\gtrsim 100\,\mathrm{EeV}$ [26]. Soon, the mass composition at the highest energies will also be tested more directly using surface detector data as well by TA [27].

The energy spectrum and mass composition of UHECRs encodes valuable information about the source environment, acceleration, and propagation of UHECRs [28, 29, 30]. For example, the suppression at the highest energies was long believed to be due to the GZK effect [31, 32], while it is now more probable to be at least partly due to the CR accelerators having reached their maximum energy around the GZK energy [33]. The ankle was previously interpreted as the pair-production dip of an all-proton composition [34], which is now ruled out by Auger composition measurements [21]. By comparing a model including source injection, propagation, and detection to data, parameters of the model can be constrained to gain information about UHECR sources. Many such fits have been conducted in the last 15 years, see e.g. [35, 36, 33, 37, 38, 39, 40, 41, 42, 43]. Usually, these phenomenological fits start with an assumption about the emission of the sources, which is often parameterized as a power-law with a cutoff function $f_{\mathrm{cut}}(E,...)$ that depends on the energy $E$ (plus additional parameters, see below):

[TABLE]

Here, the spectral index $\gamma$ and the fractions $Q_{0,A}$ of some representative elements with charge number $Z$ and mass number $A$ are free model parameters that are adjusted to best describe the measured energy spectrum and $X_{\mathrm{max}}$ distributions through a maximum-likelihood fit.

It is often assumed that UHECRs are accelerated electromagnetically, leading to a Peters cycle [44] where the cutoff energy scales linearly with $Z$ so that $E_{\mathrm{cut}}:=Z\cdot R_{\mathrm{cut}}$. The rigidity cutoff $R_{\mathrm{cut}}$ is then another free model parameter. Often, a broken-exponential cutoff function $f_{\mathrm{cut}}(E/(Z\cdot R_{\mathrm{cut}}))=\exp(1-E/(Z\cdot R_{\mathrm{cut}}))$ that sets in only when the energy $E$ becomes bigger than $Z\cdot R_{\mathrm{cut}}$ is used instead of a simple exponential cutoff for better interpretability of the value of $\gamma$. Recently, parameterizations with more variability have been tried out successfully, such as $f_{\mathrm{cut}}(E/(Z\cdot R_{\mathrm{cut}})=\text{sech}((E/(Z\cdot R_{\mathrm{cut}}))^{\Delta})$, see e.g. [43]. A value of $\Delta=2$ could for example stem from CR acceleration through magnetized turbulence [45]. An even more flexible parameterization $f_{\mathrm{cut}}(E,Z,A)\propto\exp(-E/(Z^{\alpha}A^{\beta}))$ has been investigated in [46], which could account for energy losses in the source environment or beyond-standard-model physics. It was shown that the Peters cycle ($\alpha=1$ and $\beta=0$) may not be the optimal choice, but that signatures of different cutoff functions are not significantly differentiable from influences of additional source populations with current data.

Usually, one homogeneously distributed extragalactic source population is assumed to dominate above the ankle. Directly below the ankle, the mass composition is relatively light. Models where the transition from Galactic to extragalactic sources occurs at these energies can however be excluded by the low level of anisotropy contrary to expectations from light Galactic CRs at ankle energies [47]. Thus, current models often assume a second extragalactic population which consists of protons only and dominates below the ankle. This proton component can be either fitted freely [42, 48], or be assumed to be directly related to the high-energy population and generated through photo-disintegration in the source environment [37, 49, 50] (sometimes related to as the "UFA"-model [37]). Explaining both components by a single population has the advantage that the smoothness of the energy spectrum, without any dips or abrupt breaks, is explained without finetuning of the normalizations of the two contributions. Examples of works where UHECR (and neutrino) production in specific types of source candidates is modeled and compared to data are [51, 52, 53, 54, 55, 56] for gamma-ray bursts (GRBs), [57, 58, 59, 60] for active galactic nuclei (AGNs), [61, 62, 63] for tidal disruption events (TDEs), [64] for starburst galaxies, [65] for black hole jets, and [66] for binary neutron star (BNS) mergers. In addition to the secondary component from in-source interactions in the high-energy population, a subdominant intermediate [42] to heavy [37] Galactic component is needed to explain spectrum and composition below the ankle.

Another way to describe the UHECR data is by a second extragalactic population dominating directly below the ankle that also emits a mixed composition, in addition to the high-energy extragalactic component [42]. In that case, no Galactic component is needed for $E>5\,\mathrm{EeV}$. In Fig. 1, the spectrum is shown, left for the model with a secondary proton component from in-source interactions [37], and right for the model with a second mixed extragalactic component [41].

Even though all aforementioned models put emphasis on different details, several solid conclusions about the source population dominating above the ankle have emerged as common between fits to spectrum and composition data:

• •

A not too strong source evolution: The redshift-evolution of the source population is often parameterized as $\propto(1+z)^{m}$ where $z$ is the redshift and $m$ a free model parameter. It has been shown that a strong source evolution with $m\gtrsim 5$, where the sources are predominantly far-away, leads to an overproduction of low-energy secondaries and hence an overshoot of the spectrum below the ankle [41, 42, 67]. Additionally, strong source evolutions can lead to an overproduction of $\gamma$-rays not compatible with current limits [41, 68], and can also overshoot the limits on cosmogenic neutrinos for the case of a secondary proton-producing population [69, 70]. A source evolution with $m=5$ is associated with intermediate-luminosity AGNs, and even stronger evolutions with high-luminosity AGNs [71], thus disfavoring both as the sole sources of UHECRs.

• •

A hard emission spectrum $\gamma<1$ in combination with a mixed (often Nitrogen-dominated) composition: This is necessary in order to describe the pronounced features of the spectrum, and the progressively heavier composition above the ankle in combination with the small mixing between elements indicated by the $X_{\mathrm{max}}$ distributions [12]. As visible in Fig. 1, each nuclear component only contributes to a small energy range. In those models, the instep is explained by the transition from Helium to Nitrogen and no additional flux contribution, e.g. from a local source, is necessary111Note that a scenario where the instep is generated by few foreground sources is disfavored also by the fact that the spectrum feature is consistently observed over the whole declination range covered by Auger [72].. The major point of criticism about these models is that values for the spectral index $\gamma\ll 2$ are unexpected from shock-acceleration. Note however that the spectrum given in eq. 1 relates to the emission leaving the source environment, so that the true acceleration spectrum can differ due to magnetic confinement and interactions in the source environment [37]. Also, the inferred value of $\gamma$ is strongly influenced by the shape of the cutoff function [43], systematic uncertainties [42], and the assumed source evolution (anti-correlation between $m$ and $\gamma$ [67]). Additionally, the flux suppression of low-energy particles in the EGMF can have a substantial effect on $\gamma$ [43, 39] - although extremely strong magnetic fields between the Milky Way and the first sources of around $B_{\mathrm{rms}}\approx 10-200$ nG have to be assumed to get $\gamma$ in accordance with shock acceleration. A consequence of the hard Peters cycle source emission is that the predicted rigidity of UHECRs above the ankle stays relatively constant and with only a small spread at $R\approx 4\pm 3\,\mathrm{EV}$ [73].

• •

Almost identical maximum source rigidities: In [74], it was investigated to what extend it is justified to approximate all sources to be the same - while even within one source class variations between source luminosity, size, magnetic field etc. are expected. They built a model where the maximum source rigidity varies between candidates in the population $\propto R_{\mathrm{cut}}^{-\beta_{\mathrm{pop}}}$. It was found that values $\beta_{\mathrm{pop}}\gtrsim 4$ are preferred, meaning that the population variance is surprisingly low within a factor of a few, and that UHECR sources are essentially "standard candles". This is necessary to explain the sharp features in the energy spectrum along with the small allowed mixing of the mass composition. In [75], the small allowed variance of the maximum rigidity is confirmed, while they found that the spectral index allows for a larger variation. Explaining the unexpectedly narrow maximum rigidity range is currently one of the pressing question about UHECR sources. One proposed explanation is that BNS mergers are the sources of UHECRs, for which the small variance in NS masses driving the dynamo gravitationally could explain the similarity between maximum rigidities [76].

Having established these findings, recent works have started to extend the combined fit framework, e.g. to constrain parameters of a possible EGMF [43, 73], or the contribution of ultraheavy elements [77]. An exciting extension is also the inclusion of explicit source candidates in the model whose contribution to the flux can then be determined. Examples include FR0 radio galaxies [78], specific AGNs [59, 60], catalogs of starburst galaxies, AGNs, or the nearby radio galaxy Centaurus A (hereafter Cen A) [42], or sources with different number densities following the large-scale structure (LSS) [73]. As in most of these analyses the arrival directions have been used as an additional, powerful observable, these works will be discussed in more detail below.

In the future, the combined fit of spectrum and composition will become even more constraining thanks to the inclusion of composition data from the surface detectors of Auger and TA as described above. The strong increase in statistics and the addition of composition information at the highest energies will lead to stronger constraints on the source parameters and model variations.

3 Constraints from cosmic-ray arrival directions

On top of the constraints from the UHECR spectrum and composition measurements, the arrival directions offer the possibility to learn more especially about the spatial distribution of UHECR sources. As the arrival directions are heavily influenced by the Galactic magnetic field, its effect on the arrival flux of UHECRs will be described in the following section 3.1. Afterwards, the constraints from large- and intermediate-scale anisotropies as well as the highest energy events will be analyzed in sections 3.2, 3.3, and 3.4, respectively.

3.1 Influence of the Galactic magnetic field

The Galactic magnetic field can lead to substantial deflections of UHECRs that have to be accounted for in order to draw conclusions about their sources. The typical deflection of a UHECR with rigidity $R:=E/(e\,Z)=10\,\mathrm{EV}$ is around $30^{\circ}$ as shown in Fig. 2(a). Another important effect of the GMF that has to be considered is the (de-)magnification of the flux from different directions [79, 80, 81]: CRs from some directions can reach Earth easily, while CRs e.g. from behind the Galactic center are deflected strongly and simply never reach Earth. This effect is demonstrated in Fig. 2(b) (for the UF23 GMF model suite, see below). It is visible that at $R=5\,\mathrm{EV}$ large parts of the flux from behind the Galactic center and up to the Galactic north and south are strongly demagnified and hence almost invisible to us (see also [81] for examples with other rigidities). It is important to note that Liouville’s theorem [82] is not violated by the flux (de-)magnification. An example where demagnification is important is the prediction of large-scale UHECR flux features such as the dipole [81] as discussed also below in sec. 3.2. In Fig. 2(c), the influence of different GMF models (that will be described below) on an extragalactic flux distribution is shown as a function of rigidity. In the upper plot, that extragalactic flux is solely dipolar and points in the direction of the local extragalactic matter dipole component (shown as a grey square in Fig. 4(b), see [73] for details). In general, the GMF dissolves the anisotropy and with that the dipole in the extragalactic flux - an effect that becomes stronger for smaller rigidities where deflections are larger. The GMF influence on the dipolar flux component becomes more complicated when considering also the fine structure in the extragalactic flux (expected from the local matter distribution) instead of just its dipole component as shown in the lower plot. In that case, it can easily happen that the dipole is dissolved more by the GMF at rigidities $R\approx 5\,\mathrm{EV}$ than at $R\approx 2\,\mathrm{EV}$ (which is the relevant range for UHECRs as described above), or that the dipole amplitude is even intensified by the GMF through magnification (such as for the KST24 model).

There are multiple tracers for the magnetic field strength such as Faraday rotation measures or polarized intensity of synchrotron radiation that can be used for modeling the GMF, for a review see e.g. [85]. The main challenge is that the tracers are only available integrated along the line of sight, which makes the large-scale field structure especially complicated to infer due to Earth’s position within the Milky Way disk. In the last decade, the state of the art GMF model was the JF12 model [86, 87] which consists of parametric model parts for the disk field, as well as the toroidal and poloidal fields. It can be divided into a regular field part responsible for coherent deflections and a random field part. Both are of approximately equal strength. The random field part was updated when new data from the Planck satellite became available [88].

In the last two years, great efforts have been made to improve on different aspects of GMF modeling, leading to the emergence of quite some new GMF models. Note however that to this point, only the regular part of the GMF has been studied in these updates. Thus, for now all results shown in this review are produced using the JF12 Planck-tuned random field [88] (abbreviated as "Pl" in all figures).

The UF23 suite of GMF models [13] provides for the first time an uncertainty estimate of the GMF model by exploring combinations of different parameterizations, data sets used for the tracers, and auxiliary models such as the model for the distribution of thermal electrons. This has lead to a suite of 8 models which are reasonably representative of the expected uncertainty. Additionally, newer tracer data has been taken into account which is not in agreement with the JF12 model anymore. The KST24 model [14] considers for the first time the contribution of the local bubble (a magnetized region of $\mathcal{O}(200\,\mathrm{kpc})$ in radius in which the Sun resides [89]) to the tracers during modeling. Additionally, they included the Fan region (a region of strong radio emission near the Galactic plane) as a regular Galactic-scale feature instead of masking it out as local foreground - which leads to very strong inferred field strengths in the Perseus spiral arm and is still under debate [84]. After the release of the KST model, another (preliminary) version within the UF23 model suite was produced that also takes into account the effect of the local bubble, called the UF23-locBub model [84]. It results in a slightly better fit quality than the previous 8 UF23 models. Note that even when including a more realistic shape for the bubble in variations of the UF23 model, the predicted UHECR deflections stay within the ballpark of the other UF23 models [90]. After discrepancies in the demagnification of the UF23 model suite compared to the predecessor JF12 model were found in [81] (significantly influencing the predicted large-scale UHECR anisotropy, see sec. 3.2) another model variant of the UF23 suite was introduced, the UF23-asymT model [84]. For that model, the UF23 fit was forced into a local minimum where the radial extent of the toroidal halo fields is different in the north and south, just as is the case for JF12 (and KST24). The UF23-asymT model reaches a $\chi^{2}$ within the ballpark of the other UF23 models, implying that an asymmetric halo cannot be excluded from current rotational measure data.

The deflection of a UHECR with $R=20\,\mathrm{EV}$ within all aforementioned GMF models is depicted in Fig. 2(d). In general, a relatively good level of agreement between the models is visible which is promising for UHECR source studies. The JF12 model generally predicts deflections within the range of the UF23 models. The deflections by the KST24 model are often in that region as well due to the shared antisymmetric halo, but are generally larger (as expected from Fig. 2(a)). This is especially true in the outer Galaxy due to the strong model field strength in the Fan region.

3.2 Constraints from large-scale anisotropies

The only significant anisotropy in the UHECR arrival directions is a large-scale dipole above 8 EeV [91]. It has been measured with a current significance of $6.8\sigma$ [92] as a right-ascension modulation at an energy $E>8\,\mathrm{EeV}$ by Auger. Using the full-sky data set by Auger+TA [93], its amplitude can be constrained without assumptions about the higher moments, to $6.5\%$ [94]. TA data alone is compatible with both isotropy and a dipolar distribution, with a slight preference for the latter [95]. The dipole amplitude increases with the energy as shown in Fig. 3 (left), and the amplitudes are compatible between the Auger-only and Auger+TA measurements. Note that the dipole moment is the only significant one, and that all higher multipole moments are compatible with isotropy as shown in Fig. 3 (right). The dipole direction is shown in Fig. 4(a). It points $\sim 114^{\circ}$ away from the Galactic center [94], indicating an extragalactic origin of UHECRs at these energies. The dipole direction measured by Auger-only is compatible with the full-sky Auger+TA one apart from the highest energy bin where the angular difference is $\vartheta=81^{\circ}$. This is due to intermediate-scale flux excesses visible in TA data [96] (mainly the "TA hotspot" which will be discussed further in sec. 3.3) driving the direction of the Auger+TA dipole more towards the TA field of view (FOV).

Explaining the dipole has been the aim of several publications in the last decade. Different theories have emerged:

• •

A few sources dominate the UHECR flux $>8\,\mathrm{EeV}$ and generate the dipole. That scenario is for example explored in [59, 60] where a model based on a few nearby radio galaxies in combination with a diffuse background from farther-away unresolved radio galaxies is used to explain spectrum, composition, and large-scale anisotropies. In most cases discussed in [60], the EGMF has to be rather strong with around $1\,\mathrm{nG}$ for $1\,\mathrm{Mpc}$ coherence length in order not to produce too strong anisotropies, and the flux is dominated by Virgo A (situated in the Virgo cluster) and Fornax A (see also [97]). With only few dominant sources, reproducing the subdominant quadrupole, the dipole direction energy evolution, and the spectrum features that are visible over the whole sky [72, 19] requires finetuning.

• •

An even more extreme case would be that only one source is responsible for the dipole, a scenario that was explored already long before high-statistics UHECR measurements and models of the GMF were available, see e.g. [98, 99]. It is an especially compelling hypothesis because the dominance of a single source would naturally explain the observed similarity of the sources described in sec. 2 as well as the directional uniformity of the spectral features. Using newest measurements of the spectrum and composition, in [100] a model with the nearest radio galaxy Cen A (at around $3.8\,\mathrm{Mpc}$ distance) supplying the flux above the ankle is tuned to the data. Qualitatively, the model reproduces the large-scale anisotropy for a very strong EGMF field strength of around $(20-50)\,\mathrm{nG}$ that has to be present between Cen A and Earth. If such models can also quantitatively reproduce the arrival directions over the whole energy range will be very interesting to see in the future [101].

• •

UHECR sources are numerous and at least roughly follow the large-scale structure (LSS) as explored in [102, 103, 104, 105, 91, 106, 107, 108, 109, 110, 111, 73, 81, 92]. These sources could be compact objects tracing the matter distribution, transient events occurring more often in matter-dense regions, or the acceleration of UHECRs could happen in accretion shocks present around galaxy clusters and filaments [112]. The case where UHECR sources are numerous is favored by the fact that the UHECR spectrum is smooth without any bumps hinting at contributions by individual sources. In this scenario, the dipole is generated by sources following the anisotropic local matter distribution dominated by galaxy clusters, most importantly Virgo, Great Attractor, and Coma (all in the Galactic north), as well as Perseus-Pisces. This is visible in Fig. 4(b) showing the flux at the edge of the Galaxy expected from such a model. Note that in such models, the dipole is mostly generated by (sources in) the galaxy clusters in the Galactic north, whose flux is then coherently deflected southwards to reproduce the observed dipole direction (see Fig. 4(a)). The absence of a flux excess in the Virgo direction that is sometimes regarded as peculiar (e.g. [113]) is hence naturally explained by GMF deflections. It was demonstrated in [73] that the (dark) matter distribution can be used as a bias-free estimator of the UHECR source distribution, meaning that neither an increased amount of sources in overdense nor in underdense regions is preferred (see also [102]). Especially if Virgo as an overdense cluster region is not emitting UHECRs (e.g. due to magnetic confinement [114]), the dipole cannot be well reproduced with that model. As the case of multiple sources following the LSS is a natural assumption that is frequently discussed in the literature, and predictions about anisotropies and constraints on quantities such as source density and distribution have been drawn using the newest models of the GMF, the following part of this subsection will concentrate on that scenario.

For the large-scale anisotropies to be explainable by sources following the LSS, the composition has to be mixed because a protonic composition is not in agreement with the measured dipole direction [110] and also disfavored by too-large quadrupole moments [107]. In [73], a model where the source distribution is assumed to follow the (dark) matter distribution as provided by CosmicFlows2 [115] is fit to the spectrum, composition, and dipole moments $>8\,\mathrm{EeV}$. The model source emission parameters are in agreement with the findings described in sec. 2. The resulting contribution of different element groups from different distances is shown in Fig. 5. Above 8 EeV, the flux is dominated by $\sim 53\%$ Helium and $\sim 41\%$ from the CNO group, while at $>32\,\mathrm{EeV}$, $\sim 58\%$ is CNO, $\sim 29\%$ Si-like and $\sim 12\%$ Fe-like. The effect of different horizons of different elements is visible in Fig. 5: at these energies, lighter elements like Helium have much smaller interaction lengths than heavier elements, leading to a mass ordering in the distance from which different elements contribute [116]. Note that that also implies a significantly larger dipole moment of lighter elements simply from propagation because the anisotropy is generated by the first $\lesssim 200\,\mathrm{Mpc}$ after which the universe becomes more and more homogeneous on large scales. In this model, the dipole stems mostly from primaries, moving from a dominant Helium contribution at 8 EeV to mostly CNO at 32 EeV. Secondaries, decay products from heavier primaries, typically come from further away as visible in Fig. 5 and are hence more isotropically distributed.

Another important consequence of this model is that the observed rise of the dipole amplitude with the energy is naturally explained by the shrinking propagation horizon with the energy: only around 25% of the flux comes from within the more anisotropic $100\,\mathrm{Mpc}$ above 8 EeV, which increases to 65% above 32 EeV as visible in Fig. 5. As explained in sec. 2, the rigidity stays almost constant with the energy as the mass composition gets heavier, so that the suppression of the dipole for lower rigidities by the GMF as shown in Fig. 2(c) is only a subdominant effect in this model. The evolution of the dipole and quadrupole moments predicted by the LSS-based model from [73] for different GMF models was explored in [81] and can be seen in Fig. 6 (upper row). It is visible that for the continuous model (in the limit of infinite source number density), the dipole amplitude is consistently smaller for the UF23 models compared to the predecessor JF12 model. This is because the Virgo/Great Attractor overdensity region is demagnified (Fig. 2(b)) by all UF23 models compared to JF12 [117, 81], as visible in Fig. 2(c). The underlying reason for this is the symmetric toroidal halo field of the UF23 model variants [84] - when using the UF23-asymT model with forced asymmetric halo, the values are indeed more similar to the ones predicted by JF12 as visible in Fig. 6. The KST24 model also has an asymptotic halo and magnifies the Virgo/Great Attractor region (see Fig. 2(c)), leading to a larger dipole component.

When instead of the continuous model, the more appropriate case of a finite source number density is investigated, the dipole and quadrupole moments are influenced by cosmic variance: they vary depending on the exact source locations, which is indicated by the $1\sigma$ error bars in Fig. 6. The smaller the source density, the larger the contribution by less and less individual local sources becomes, driving up the anisotropy. For a source density of around $n_{s}=\mathcal{O}(10^{-4}\,\mathrm{Mpc}^{-3})$, the best agreement with dipole and quadrupole moments is reached for the UF23 models, as larger densities undershoot the dipole and smaller ones overshoot both dipole and quadrupole. For the older JF12 model, the KST24 model, and the UF23-asymT model, only larger source densities $n_{s}\gtrsim 10^{-3.5}\,\mathrm{Mpc}^{-3}$ are compatible with both quadrupole and dipole amplitudes. This indicates how sensitive the inferred source density is to properties of the GMF model - and that the number has to be treated cautiously, see also the discussion below. Note that these results based on the CosmicFlows 2 LSS model are compatible with studies based instead on the 2MRS galaxy catalog [111, 92]. The obtained limits on the source density allow for conclusion on the sources of UHECRs. Sparse source types like blazars ($n_{s}\sim 10^{-7}\,\mathrm{Mpc}^{-3}$ [118]), high-luminosity AGNs ($n_{s}\sim 10^{-6}\,\mathrm{Mpc}^{-3}$ [119]), and even starburst galaxies ($n_{s}\sim 10^{-5}\,\mathrm{Mpc}^{-3}$ [119]) may be too sparse to be the sole sources of UHECRs due to too strong anisotropies if no severe smoothing by the EGMF is considered [73, 66]. More common source types like starforming galaxies, low-luminosity AGNs, or transients with high rates like non-jetted TDEs [120] or BNS merger [76] are instead preferred.

The predicted dipole directions of the CosmicFlows-based model are shown in Fig. 7. Note that the dipole directions predicted using the 2MRS catalog [92, 111] agree with the CosmicFlows predictions shown in Fig. 7 within $\mathcal{O}(10^{\circ})$ for the JF12 model. Predictions for other GMF models based on the 2MRS catalog are not yet available. The dipole directions predicted using the full sky and using only the limited Auger FOV (both depicted in Fig. 7) can differ by up to $\vartheta=30^{\circ}$ at lower energies, depending on the GMF model. For the highest energy bin $E>32\,\mathrm{EeV}$, the difference is $\vartheta\lesssim 10^{\circ}$ for all models. All GMF model predictions are roughly compatible with the measured dipole directions by Auger (even though no model fits perfectly for all energies). However, none of the models can reproduce the full-sky measured dipole direction by Auger+TA222Note that the measured $\vartheta=81^{\circ}$ angular difference between the full-sky and Auger FOV dipole directions in the $E>32\,\mathrm{EeV}$ bin is also not reproduced when considering cosmic variance. On average, it is$\sim 20^{\circ}$ for $n_{s}=10^{-3}\,\mathrm{Mpc}^{-3}$ and $\lesssim 30^{\circ}$ for $n_{s}=10^{-4}\,\mathrm{Mpc}^{-3}$, independently of the GMF model. in the highest energy bin $>32\,\mathrm{EeV}$, that differs substantially from the Auger FOV one. This is because the intermediate-scale anisotropies measured by TA at energies $\gtrsim 40\,\mathrm{EeV}$ (see sec. 3.3) pulling the dipole direction closer to the TA FOV are not reproduced by the LSS with any GMF model. Thus, if these intermediate scale anisotropies turn out to be statistically significant, they are more likely to originate from local or transient sources instead of the galaxy clusters represented in the LSS, see also [121].

The difference between predictions using different random field coherence lengths is $\mathcal{O}(15^{\circ})$ at lower energies and decreases to $\mathcal{O}(5^{\circ})$ for $E>32\,\mathrm{EeV}$ (not shown, see [81]). The predicted dipole direction clearly also depends on cosmic variance, and variations increase fast with decreasing source density. In Fig. 7, the $1\sigma$ uncertainty for the UF23-base model at $n_{s}=10^{-3}\,\mathrm{Mpc}^{-3}$ is shown exemplarily. Already at that source density, the uncertainty from cosmic variance is larger than the differences between different GMF models.

Next to the uncertainty from cosmic variance which is an inherent uncertainty that is not avoidable until the sources of UHECRs are identified individually, there are evidently other uncertainties involved in models such as the ones described above. The uncertainty on the coherent GMF model is manageable thanks to the recent development of several new models [81]. Remarkably, even GMF models by different authors agree on both the dipole direction and amplitude predictions within the uncertainties from cosmic variance for source density $\lesssim 10^{-2}\,\mathrm{Mpc}^{-3}$. The EGMF, however, could have an unknown effect, both on the inferred source density and also on other model parameters as investigated in [73, 39, 43, 111]. But, to estimate that effect reliably, more knowledge about the structure and strength of the EGMF [15, 122, 123], that of nearby galaxy clusters [114], and the magnetic field of the local group are necessary. When considering also structured EGMFs instead of simple turbulent approximations, simulations and model fits become extremely computationally expensive.

Another important uncertainty stems from the source model - while it seems that conclusions reached with models based on LSS approximations from dark-matter distributions such as CosmicFlows 2 [106, 109, 110, 73, 81] agree in general with the conclusions drawn from models based on galaxy catalogs [107, 92, 111], it was also shown that the predicted large-scale anisotropies are extremely sensitive to the interplay of GMF model demagnification and source distribution. Updated predictions using CosmicFlows 4 [124] as well as explicit galaxy catalogs with newer GMF models are needed for better estimation of the impact of the source distribution. Note that also a homogeneous source distribution can lead to the observed large-scale anisotropies for sufficiently small source densities $n_{s}\lesssim 10^{-4}\,\mathrm{Mpc}^{-3}$ [104, 108, 125, 111, 73], but for that case the dipole direction is clearly non-informative.

Finally, the mass composition can also have significant impact on the predicted large-scale anisotropies. While variations of the hadronic interaction model or source emission fit parameters seem to have minor impact [73, 111], significant shifts of the $X_{\mathrm{max}}$ scale towards heavier elements could lead to a decreased level of anisotropy [126]. The measurement of mass estimators for individual events at the highest energies [22] also offers great chances to learn more about UHECR sources and the GMF. E.g. the analysis proposed in [127] where the dipole component is calculated for a heavy and a light event selection in each energy bin could test the mass-dependent effects described above: the anticorrelation between dipole moment and mass expected from propagation, and the rigidity-dependent suppression by the GMF as shown in Fig. 2(c). Especially a significant split of the dipole direction between heavy and light selections could be extremely informative about GMF models and possibly even the influence of local sources.

3.3 Constraints from intermediate-scale anisotropies

While at energies $\sim 8\,\mathrm{EeV}$ the UHECR sky is dominated by the large-scale dipole signal, at higher energies $\gtrsim 32\,\mathrm{EeV}$ anisotropies at smaller scales start arising. This is as expected from the shrinking propagation horizon as visible in Fig. 5, indicating a larger relative flux contribution of local sources. The full-sky flux at energies $\gtrsim 40\,\mathrm{EeV}$ as measured by Auger and TA combined is shown in Fig. 8 (left). The Auger data has been scanned for anisotropies over different angular scales and energy thresholds in every available direction, revealing no significant anisotropies beyond a $p$-value of $2\%$ [128, 129]. The maximum significance is reached at $E>38\,\mathrm{EeV}$ for an angular scale of $27^{\circ}$. It is close to the direction of the nearby radio galaxy Cen A, which will be discussed further below. The autocorrelation and correlation with the Galactic and supergalactic plane all yield $p$-values $>10\%$. In the TA data, a scan has revealed an overdensity at $E>57\,\mathrm{EeV}$ which is often referred to as the "TA hotspot" [96]. Its significance grows over time and is now at $2.9\sigma$ post-trial on an angular scale of $25^{\circ}$ [130]. Several theories for its origin have been discussed in the literature, mostly nearby candidates close to the hotspot direction [121] like the starburst galaxy M82, or more complicated proposals like CRs originating in Virgo and then being deflected along filaments [131] (note however that a pure-proton composition was used in [131]). Another excess at slightly lower energies $\sim 25\,\mathrm{EeV}$ was later identified in a manual search at $20^{\circ}$ angular windows with a local significance of $\sim 3.9\sigma$. Often, this overdensity is referred to as the "Perseus-Pisces" (PP) excess because it is close to the supercluster direction (compare to Fig. 4(b)). Note that the quoted post-trial significance of $3.2\sigma$ [130] refers to how often such an excess appears close to PP instead of quantifying the overdensity on a statistical basis accounting for the scan in direction and energy that was indirectly performed. Note also that in the studies described above based on a LSS model including PP as a dominant source in the extragalactic flux, its flux excess is dissolved by the GMF and not observable at Earth [73].

The combined full-sky significance map from Auger and TA is shown in Fig. 8 (right) with the names of the three described intermediate-scale overdensities indicated. Interestingly, both the TA hotspot and the PP-excess are also within the FOV of Auger, albeit at the border and only observed in inclined events. However, due to the larger size of Auger, its exposure in the two excess regions is comparable to the one of TA. In Auger data, no overdensity in neither the hotspot nor the PP direction is observed [128, 132]. This also disfavors the claim [133] that flux excesses in the northern hemisphere only visible for TA could explain the discrepancy in the flux seen by Auger and TA, a hypothesis investigated e.g. in [134, 135]. Note that taking into account also possible energy scale mismatches, a mere statistical overfluctuation seen by TA, or an underfluctuation seen by Auger in the excess regions cannot be excluded at this point [132].

Correlation with Cen A:

The correlation with the nearby radio galaxy Cen A is currently observed with $4\sigma$ level post-trial significance for $E>38\,\mathrm{EeV}$ [128] at an angular scale of $27^{\circ}$. Cen A has been suspected as a source of UHECRs for years [136, 137, 138, 139], even long before the overdensity had emerged, due to its proximity and powerful radio jets. Acceleration of UHECRs up to the highest energies is possible for Cen A according to simulations and theoretical modeling, e.g. [140, 141].

In [42], a model was simultaneously fit to the Auger energy spectrum, mass composition, and arrival directions $>16\,\mathrm{EeV}$, based on the assumption that Cen A (or other galaxy catalogs, see below) is a source of UHECRs and contributes part of the flux above the ankle. Additionally, another part of the model flux originates from homogeneous background sources, as in the models introduced above in sec. 2. By taking into account propagation effects and a rigidity-dependent blurring (but no coherent deflections), the level of anisotropy rises naturally with the energy due to the increased relative contribution of the nearby source Cen A with a best-fit of around $3\%$ at $40\,\mathrm{EeV}$. The model including Cen A describes mostly the arrival direction distribution significantly better than the model based on only homogeneous sources, especially around $10^{19.3}\,\mathrm{eV}\approx 20\,\mathrm{EeV}$ and $\sim 10^{19.7}\,\mathrm{eV}\approx 50\,\mathrm{EeV}$. These two peaks in the correlation could hint at different mass groups being accelerated to the same rigidity [142]. In the model discussed in [42], the dominant mass groups at the respective energies are Helium and CNO with an appropriate factor 2-3 between their charge numbers.

In [132], it was recently confirmed that the Cen A excess extends to smaller energies of $\sim 20\,\mathrm{EeV}$. The position of the excess remains constant with the energy, which could indicate that the CR rigidity remains constant or that no strong coherent deflections occur due to a subdominant regular magnetic field in that direction. Note that the angular window was kept fixed in the analysis [132], so no conclusions regarding the GMF can be drawn from that. The arrival directions of UHECRs from Cen A that are expected for different GMF models are shown in Fig. 9. It is visible that all GMF models predict that CRs are deflected towards the Galactic plane, which is also the direction where the maximum overdensity occurs, see the cross marker in Fig. 9. Overall, coherent displacement remains rather small, for all GMF models apart from the KST24 model. Comparing these simulations visually to Fig. 8, it does seem reasonable that Cen A could be the origin of the flux excess if at least part of its flux consists of events with $R\gg 4\,\mathrm{EV}$, possibly in combination with some additional blurring from a turbulent EGMF that is not considered in the simulations shown in Fig. 9, see also [143] and [101].

Correlation with active galactic nuclei:

Not just the nearest radio galaxy Cen A, but also active galactic nuclei (AGNs) in general were long suspected to be sources of UHECRs due to their powerful jets and magnetized lobes, both of which could be good candidates for acceleration to the highest energies [7, 144]. Different theories for acceleration in the core, lobes, and jets of AGNs (including also two-stage and re-acceleration processes) have been investigated, see e.g. [145, 146, 147, 148, 149]. In the case that UHECRs are accelerated along the jet, blazars where the jet points into the direction of Earth are prime candidates for UHECR acceleration [150, 57]. Note that the detection of a correlation of the UHECR arríval directions with AGNs might be hindered if time delays e.g. due to the EGMF become larger than typical AGN duty cycles [151].

Early Auger data $>60\,\mathrm{EeV}$ seemed to correlate with AGNs on angular scales $\lesssim 6^{\circ}$ [152]. While that correlation turned out to be a statistical fluctuation, current tests for correlation between the Auger+TA arrival directions and AGNs have reached a significance of $3.3\sigma$. Three parameters were optimized, the energy threshold $E>38\,\mathrm{EeV}$, an anisotropic fraction $4.8\%$ and a Fisher search radius around the sources of $15.4^{\circ}$ [153]. Additionally, a model of jetted AGNs ($\gamma$-AGNs) has been tested, where the UHECR luminosity is assumed to be scaled with the $\gamma$-ray flux measured by Fermi. The model reaches a significance of $\sim 3.8\sigma$ for a similar angular scale and energy threshold and an anisotropic fraction of $8.8\%$. In the subsequent analysis including an energy-dependent modeling of the catalog contribution in a combined fit to spectrum, composition, and arrival directions measured by Auger, the $\gamma$-AGN model was however disfavored due to the strong dominance of the far-away blazar Markarian 421 [42]. In [154] it was suggested to correct the UHECR flux weighting with the observed $\gamma$-ray flux, as that implicitly assumes that UHECRs experience the same beaming along the jet as $\gamma$-rays - which is unexpected due to magnetic field decollimation. Instead, for example the intrinsic unbeamed $\gamma$-ray flux or the jet power could be used as a proxy. This is also in accordance with the finding from [59, 60] that nearby radio galaxies could explain the measured Auger data including the dipole when the UHECR flux is more carefully modeled and not simply scaled with the $\gamma$-ray flux (see sec. 3.2), and the observation that $\gamma$-ray bright sources like blazars alone are too sparse and hence produce too much anisotropy [155, 73]. Instead of considering only the brightest and most powerful AGNs that lead to too large anisotropies, lately also models considering weaker more prevalent classes of AGNs like FR0 [156, 78] or Seyfert galaxies [157] have been investigated, even though it remains unclear if they can produce sufficiently high energies. Those weaker source candidates also have the advantage that their source evolution is weaker [71] and hence does not lead to the overproduction of secondaries overshooting the spectrum below the ankle as is the case for powerful AGNs [41, 42], see sec. 2.

Correlation with Starburst galaxies:

Next to AGNs, also a catalog of starburst galaxies (SBGs) has been tested for correlations with arrival direction data. These are galaxies detectable by high infrared luminosities as they undergo extreme starformation activity, driving powerful magnetized winds. Comparing the Auger+TA arrival direction data to a model where part of the UHECR flux comes from these SBGs, currently a significance of $4.2\sigma$ is reached for an energy threshold $E>38\,\mathrm{EeV}$, an anisotropic fraction $\alpha=10.6\%$ and a search radius of $17.6^{\circ}$ [153]. Using again the model including propagation and rigidity-dependent blurring for a combined fit of the Auger energy spectrum, composition, and arrival directions, the SBG model is favored over a homogeneous model at $4.5\sigma$ significance (including experimental systematic uncertainties) [42]. Especially the excess in the Cen A region is well described by the SBG model. In that model, the flux excess comes from NGC4945 instead of Cen A, an SBG close in direction and at a similar distance $\sim 4\,\mathrm{Mpc}$.

Despite its quite large significance, the observed correlation with SBGs is difficult to interpret because only symmetric search radii or a rigidity-dependent blurring was taken into account in the analyses [153, 158, 42], but no coherent displacements as expected by the regular part of the GMF. In [117, 159], it was shown that it is possible but not straightforward to reproduce the SBG correlation if UHECR sources simply follow the 2MRS catalog and GMF deflections are included. The correlation is however easier to reproduce in models with coherent deflections if UHECR sources are galaxies with high starformation rate [117]. In [160, 161], it was also shown that it is possible to reproduce the significance, the best-fit search radius, and the anisotropic fraction inferred by Auger in simulations based on the SBG catalog even when including a coherent GMF model. In that case, the true signal fraction from the SBG catalog has to be higher $\sim(10-30)\%$ (depending on the GMF model) than the inferred anisotropic fraction of $\sim 10\%$ because heavier particles are deflected too much and do not fall within the search radius anymore [161].

Naturally, the observed correlation between the UHECR data and starburst galaxies has sparked debate if SBGs could actually be the sources of UHECRs [162, 163, 164]. It is still somewhat unclear if and how CRs could be accelerated to ultra-high energies in SBGs, especially because strong radiation losses are expected [165, 166, 167, 168, 64]. Alternatively, UHECRs could also be (pre-)accelerated not by the SBGs themselves, but by transient events like supernovae or GRBs which are very numerous in the central region of SBGs [169]. The hypothesis that transient events occuring proportionally to the SFR in any galaxy are the sources of UHECRs is investigated in [170]. By comparing models with different transient burst rates to the flux excesses seen in the Cen A region and the TA hotspot region, they conclude that long GRBs are favored to explain the measured data. Note however that that conclusion may change once coherent magnetic field deflections and event-by-event data is taken into account [171]. Another possibility to explain the observed correlation of UHECRs with SBGs without SBGs being the actual sources of UHECRs is explored in [172, 173]. In that model, UHECRs are accelerated by Cen A in an earlier powerful outburst episode and then merely scattered by magnetic fields of the nearby galaxies forming the "Council of Giants" - which includes many of the SBGs used in the correlation analysis described above [158, 129].

Deflection patterns and multiplets:

While many of the described analyses hint at nearby source candidates like Cen A or SBGs being sources of UHECRs, no firm conclusions can be drawn yet. Even if a correlation significance surpasses $5\sigma$, it could still be that the overdensity is actually caused by other sources. For the Cen A excess for example, other source candidates next to Cen A, such as the nearby starburst galaxy NGC4945 (see above), or further away sources like the Virgo cluster could also explain the observed anisotropies [110, 111, 117].

One way to disentangle the different sources in the Cen A excess region would be to search for energy-ordered events (multiplets) as expected from magnetic field deflections. In a blind search, no significant multiplets were found in the Auger FOV [174, 175]. Also, no elongated patterns expected from coherent deflections were found in a principal-axis analyses. In a targeted search at Cen A’s position, a 9-plet was found with a $p$-value of 6% [175]. In another preliminary scan on Auger open data using a dedicated likelihood including parameters of the magnetic field [176], a multiplet is identified with a post-trial significance of $3.3\sigma$, including also one of the highest-energy events measured by Auger at $165\,\mathrm{EeV}$. Note that the multiplet itself is in the Cen A excess region, yet - considering coherent magnetic field deflections - the most probable source is identified to be the Sombrero galaxy, an AGN at around $9\,\mathrm{Mpc}$ distance. Note that the likelihood to find multiplets will in general always be higher in flux excess regions (such as the Cen A region). This effect is not penalized for in [176]. They also identify a multiplet pointing towards Cen A, albeit with a lower significance. Another targeted search conducted by TA has found a $4.1\sigma$ indication for energy-ordered multiplets along the supergalactic plane [177] that could be indicative of UHECRs diffusing from sources along that plane and hence roughly the LSS.

In the future, it will be very insightful to search for rigidity-ordered multiplets instead of just energy-ordered ones using the event-level mass estimators being developed by Auger [22] and TA [27]. This is especially important because the UHECR mass composition changes so quickly with energy. A rigidity-ordered multiplet pointing unambiguously into the direction of a source like Cen A would be a clear indication of it being an actual UHECR source, and also the first direct observation of coherent magnetic field deflections. On the other hand, the reconstruction of the magnetic field indirectly from UHECR deflections might also become possible [178], especially using machine-learning methods (see e.g. [179, 180]).

3.4 Constraints from the highest-energy events

At even higher energies than where the intermediate-scale anisotropies have started to arise, the event number decreases quickly. But, due to the shrinking propagation horizon with the energy, these events have to originate from relatively nearby sources so that correlation studies with nearby matter or auto-correlation may still be very informative. Before the start of Auger and TA, when the event statistic was still very limited, the arrival directions of UHECRs $>40\,\mathrm{EeV}$ seemed to exhibit clustering in doublets and triplets [181, 182], which was utilized to constrain properties of UHECR sources like the source number density, see e.g. [183, 184, 185, 186, 187]. Using an early Auger data set from 2011, a limit of $n_{s}>6\times 10^{-6}\,\mathrm{Mpc}^{-3}$ was set on the source number density $>70\,\mathrm{EeV}$ - but that limit only applies if magnetic field deflections are $<30^{\circ}$. As visible in Fig. 2(a), deflections can be a lot larger for a realistic mass composition heavier than nitrogen at $70\,\mathrm{EeV}$. No comparable analyses taking into account heavier mixed composition and updated GMF models are yet available.

Recently, TA has measured an event with an exceptionally large energy of $E=244\pm 29\,(\mathrm{stat})^{+51}_{-76}\,(\mathrm{syst})\,\mathrm{EeV}$ referred to as the "Amaterasu" event [188]. By backtracking its arrival direction using different models of the GMF and assuming that the particle was originally an iron nucleus, it was found that Amaterasu most likely originated from a direction within the local void, and cannot be associated with any known local powerful galaxy or cluster [188, 189, 83, 190]. This is shown in Fig. 10 for different models of the GMF and associated uncertainties. Due to its large energy, the source has to be close [191], within $\sim(8-50)\,\mathrm{Mpc}$ depending on the exact event energy and assuming an iron nucleus at the source [189]. For all TA events $>100\,\mathrm{EeV}$ to be associated with the local matter distribution, either a very strong EGMF has to be present, or the composition has to be extremely heavy [26]. For Auger, a statistical analysis was performed to test the association between the events $>100\,\mathrm{EeV}$ and different source catalogs [192], as displayed in Fig. 10 (right). Not even a small contribution from any of the catalogs fits significantly better than isotropy. A contribution $\gtrsim 60\%$ from the Lunardini catalog of starburst galaxies, the Swift-BAT AGN catalog, and the Fermi-Lat catalog of jetted AGNs is even excluded at $5\sigma$. Note that deflections by the EGMF are currently neglected in the analysis.

This non-association of the highest energy events with local source candidates or the LSS has sparked theories about the origin of Amaterasu and the highest energy events in general: next to a strong EGMF [100] or more exotic theories like beyond-standard-model physics or superheavy dark matter [193, 194, 195], it could simply be that the sources indeed lie in voids and are just not special or powerful enough to be a prominent source candidate. In that case, UHECRs could originate from transient events happening in these ordinary galaxies [189], such as TDEs, young magnetars, or BNS mergers [76]. Another possibility is that UHECR generation anti-correlates with the electromagnetic source power, so that less prominent sources like the Sombrero galaxy (see sec. 3.3) are responsible for UHECR production [196]. For these lower-power sources, lower photon densities are expected, which may allow UHECRs to escape without heavy losses and hence without producing a strong $\gamma$-ray flux, see also [155].

Another possibility is that the highest energy events are in fact heavier than iron [197] and produced for example in r-processes happening in BNS mergers or collapsars [77, 76, 198]. In that case, the propagation distance is still $\mathcal{O}(100\,\mathrm{Mpc})$ at $300\,\mathrm{EeV}$ [77, 199], so that the sources could be further away. Additionally, deflections in the GMF would be larger for higher charge numbers. For the Amaterasu event, this means that sources on the supergalactic plane fall within the possible backtracked directions [77].

4 Conclusions and Outlook

Since the detection of UHECRs, their origin has remained a mystery. Thanks to the large increase in statistics of data by the Pierre Auger Observatory and Telescope Array, however, several characteristics of UHECR sources have been narrowed down over the last years. Quite a lot of the findings were unexpected, such as the intermediate to heavy mass composition at the highest energies. By adapting models to best describe the measured energy spectrum, shower depth distributions, and lately also the arrival directions, the following characteristics of UHECR sources have been identified as most likely in the energy range above the ankle ($\sim(8-32)\,\mathrm{EeV})$:

• •

They are extragalactic and evolve with a source evolution weaker than $\propto(1+z)^{5}$ [33, 41, 42, 67] (see sec. 2).

• •

They emit (possibly after in-source interactions or magnetic confinement [37]) a relatively hard spectrum and mixed composition [33, 41, 42] where the maximum rigidity does not differ much between sources [74, 75] (see sec. 2).

• •

They at least roughly follow the large-scale structure [102, 103, 104, 105, 106, 107, 109, 110, 111, 73, 81, 92] and have a relatively large density $n_{s}\gtrsim 10^{-4}\,\mathrm{Mpc}^{-3}$ (see sec. 3.2) [111, 73, 81, 92].

Especially the last two points have sparked debates about which kind of source class could be on one hand so common, and on the other hand so identical [76]. Also, further modeling is needed to gain information about the exact relevant energy range where the above listed findings apply, i.e. if the same source class could be responsible for the dipole, intermediate-scale anisotropies, and the highest energy events. The sources of the smaller-scale anisotropies are still heavily under debate, e.g. the source behind the flux excess in the Cen A region seen by Auger. The overdensity could be caused for example by Cen A itself [138, 139, 173], by the nearby starburst galaxy NGC4945 [158, 42], or by deflected events from the Virgo cluster [110, 117] or the Sombrero galaxy [176, 196] (see sec. 3.3). Another pressing question has emerged regarding the highest energy events, which cannot easily be associated with nearby prominent source candidates or the extragalactic matter distribution [188, 189, 83, 84]. If these events are produced by the same source class as dominant above the ankle, the sources are likely of transient nature - such as binary neutron star mergers [76] - and the mass composition at the highest energies is likely heavier than iron [197, 77] (see sec. 3.4).

Some of these questions could be answered in the near future by using mass estimators from the Auger and TA surface detectors [22, 27]. Combining event-by-event charge information at the highest energies with energy and arrival direction measurements will provide powerful constraints on UHECR sources. This could be done for example by extending combined fits of models to multiple observables such as [42], by investigating the difference between heavy and light components in large- and intermediate-scale anisotropies [127, 200], or by searching for rigidity orderings near source candidates. Additionally, multi-messenger analyses can be a powerful probe of UHECR sources. First hints about UHECR sources from neutrinos are for example the observation of neutrinos from the nearby starburst galaxy NGC1068 [201], a neutrino associated with a tidal disruption event [202], or the non-observation of neutrinos from gamma-ray bursts [203].

When comparing the findings listed in this review to the open questions raised in [6], it can be seen that several of those questions have been at least partially answered in the last six years. Yet, the next generation of UHECR experiments such as GCOS [204] or POEMMA [205] is needed to provide a further increase in statistics and precision of mass indicators, especially at the highest energies, to answer remaining questions. Additionally, further improvements in modeling have to be made, e.g. by even better models of the Galactic and extragalactic magnetic field.

Acknowledgements:

I thank Jonathan Biteau, Glennys Farrar, Foteini Oikonomou, Paul Simeon, Michael Unger, and Walter Winter - as well as all my colleagues from the Pierre Auger Collaboration - for fruitful discussions, and Silvia Mollerach and Armando di Matteo for comments on this manuscript. I also thank Alexander Korochkin and Michael Unger for providing their GMF models.