# Closing the Net on Transient Sources of Ultrahigh-energy Cosmic Rays

Sullivan Marafico¹ ID, Jonathan Biteau¹,² ID, Antonio Condorelli¹ ID, Olivier Deligny¹ ID, and Johan Bregeon³ ID

¹ Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France; biteau@in2p3.fr, antonio.condorelli@ijclab.in2p3.fr, olivier.deligny@ijclab.in2p3.fr

² Institut universitaire de France (IUF), France

³ CNRS-IN2P3, Laboratoire de Physique Subatomique et de Cosmologie (LPSC), Grenoble, France

Received 2023 December 27; revised 2024 June 15; accepted 2024 June 18; published 2024 August 21

## Abstract

The arrival directions of ultrahigh-energy cosmic rays (UHECRs) observed above $4 \times 10^{19}$ eV provide evidence of localized excesses that are key to identifying their sources. We leverage the 3D matter distribution from optical and infrared surveys as a density model of UHECR sources, which are considered to be transient. Agreement of the sky model with UHECR data imposes constraints on both the emission rate per unit matter and the time spread induced by encountered turbulent magnetic fields. Based on radio measurements of cosmic magnetism, we identify the Local Sheet as the magnetized structure responsible for the kiloyear duration of UHECR bursts for an observer on Earth and find that the turbulence amplitude must be within 0.5–20 nG for a coherence length of 10 kpc. At the same time, the burst-rate density must be above $50 \text{ Gpc}^{-3} \text{ yr}^{-1}$ for Local Sheet galaxies to reproduce the UHECR excesses and below $5000 \text{ Gpc}^{-3} \text{ yr}^{-1}$ ($30,000 \text{ Gpc}^{-3} \text{ yr}^{-1}$) for the Milky Way (Local Group galaxies) not to outshine other galaxies. For the transient emissions of protons and nuclei to match the energy spectra of UHECRs, the kinetic energy of the outflows responsible for UHECR acceleration must be below $4 \times 10^{54}$ erg and above $5 \times 10^{50}$ erg ($2 \times 10^{49}$ erg) if we consider the Milky Way (or not). The only stellar-sized transients that satisfy both Hillas' and our criteria are long-duration gamma-ray bursts.

*Unified Astronomy Thesaurus concepts: Cosmic ray astronomy (324); Large-scale structure of the universe (902); Magnetic fields (994); Transient sources (1851); Ultra-high-energy cosmic radiation (1733)*

*Materials only available in the online version of record: figure set*

## 1. Introduction

The astrophysical objects that accelerate charged particles to energies above 1 EeV ($\equiv 10^{18}$ eV) have long been elusive. Yet, the mechanisms by which these nuclei are emitted and the abundances of nuclear elements in the source environments are considerably constrained by the energy spectrum and mass composition of the ultrahigh-energy cosmic rays (UHECRs; Aab et al. 2014a; Abreu et al. 2021). Data collected over the past two decades show a gradual increase with energy in the mean logarithmic mass number of UHECRs, $\langle \ln A \rangle$ (Aab et al. 2014a; Watson 2022). These observations are consistent with the fundamental expectation that electromagnetic processes accelerate charged particles of atomic mass A to a maximum energy $E_{Z_A}^{\text{max}}$ proportional, at least approximately, to their electric charge $Z_A$ (Peters 1961). The intensity of individual nuclear components, ranging from fully ionized He to Fe, is thus expected to drop off at the same magnetic rigidity at the sources, which is estimated to $E_{Z_A}^{\text{max}}/Z_A e \simeq 1–5$ EV (see, e.g., Aloisio et al. 2014; Taylor et al. 2015; Aab et al. 2017a).

The low variance observed for the mass estimators, $\sigma^2(\ln A)$, imposes a nearly mono-elemental composition at all energies above 5 EeV. The little mixture above this so-called ankle feature of the energy spectrum is reflected in a very hard spectral index for each of the nuclear species: $-1 \lesssim \gamma \lesssim 1$, for an emission spectrum that follows $dN/dE \propto E^{-\gamma}$ with a cutoff at the highest energies. Whatever the exact value of the spectral index, which depends on the systematic uncertainties affecting the various models

(Aab et al. 2017a), the hardness of the nuclear spectra differs from expectation of diffusive shock acceleration models, $\gamma \gtrsim 2$ (see Sironi et al. 2015 for a review). Acceleration by back-and-forth bounces on approaching magnetic fields, described by the first-order Fermi mechanism, can however be reconciled with observations if the shaping of ejected spectra takes into account interactions within the sources. These interactions tend to favor the escape of high-energy nuclei relative to lower-energy nuclei (Globus et al. 2015; Unger et al. 2015; Biehl et al. 2018a; Fang & Murase 2018; Supanitsky et al. 2018). A generic prediction of such in-source interaction scenarios is a proton spectrum softer than the nuclear spectra, due to the escape from the magnetized interaction zone of neutrons that subsequently decay in protons. The requirement of a soft proton spectrum and hard nuclear spectra escaping from the sources has been shown to be compatible with the UHECR data, in particular with the proton spectrum observed at energies below the ankle (Luce et al. 2022).

To supply the observed energy density, the production rate density of UHECRs is constrained to be $\mathcal{L} \simeq 10^{54}$ erg $\text{Gpc}^{-3} \text{yr}^{-1}$ above 0.6 EeV (Luce et al. 2022). The value inferred for $\mathcal{L}$ can be used to identify the few classes of sources that radiate similarly in the electromagnetic band. Several candidates can indeed meet the energetics constraints, from sources that undergo persistent or frequent episodes of electromagnetic emission with rather low luminosity to cataclysmic high-luminosity events that occur only rarely (see Alves Batista et al. 2019 for a review). The degeneracy between source density and luminosity is addressed in this work using evidence of localized UHECR excesses over the celestial sphere.

Sky maps of arrival directions support from an observational perspective that UHECRs are predominantly of extragalactic

<figure><img src="image_1.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
origin. On the one hand, the contrast of the flux distribution on the sphere above 8 EeV is suggestive of the pattern expected from the large-scale structure of matter within a few hundred megaparsecs (Aab et al. 2017b). On the other hand, the analysis of the arrival-direction data above 40 EeV together with spectral and composition data leads to the inference of a flux contribution of ≈10%–20% from the ≈40 brightest star-forming galaxies within a radius of ≈100 Mpc (Abdul Halim et al. 2024). This is consistent with the correlation of the flux pattern of these galaxies and the UHECR arrival directions alone, as previously reported by the Pierre Auger Collaboration in Aab et al. (2018) and Abreu et al. (2022). The distance scales sampled with increasing energy are in line with the reduction of the UHECR horizon, primarily due to interaction with the cosmic microwave background (CMB), also known as the Greisen–Zatsepin–Kuzmin (GZK) horizon (Greisen 1966; Zatsepin & Kuz'min 1966).

But sky maps of arrival directions contain a potentially even greater wealth of information that may break the degeneracy between frequent low-luminosity and rare high-luminosity source populations meeting the energetics requirements. For a given burst rate, the average number of sources visible at each instant depends on the effective duration of the bursts, and therefore on the magnetic wandering experienced by the UHECRs during propagation. As a result, the average number of visible galaxies hosting transient sources decreases with energy much faster than in a steady scenario. Steady and transient scenarios should therefore give rise to different sky maps. Signatures of transient scenarios have been used in the past to infer source properties and constrain propagation by supposing the detection of UHECR clustering on small angular scales (see the seminal works of Lemoine et al. 1997; Sigl et al. 1997), to constrain the local burst-rate density assuming sources distributed uniformly in space and emitting protons only (see Murase & Takami 2009), and to interpret potential features of the all-particle energy spectrum (see Miralda-Escude & Waxman 1996).

The aim of this paper is to complement the existing constraints on the production rate density of UHECRs by making use of evidence of anisotropy on intermediate angular scales to bracket the luminosity and rate density of bursting sources. The common thread of this work is that each galaxy of the flux-limited sample of Biteau (2021) can harbor transient episodes of UHECR production at a given rate. In Section 2, we describe the source-emission model that incorporates the transient nature of the UHECR production from each galaxy, the resulting intensities on Earth, the time delays and angular deflections undergone by particles during their propagation in the intervening magnetic fields, and the magnetic effects limiting the number of observable sources. In Section 3, we infer the best rate density to reproduce the main features of the sky map reported by the Pierre Auger and Telescope Array collaborations above ≈40 EeV (Abdul Halim et al. 2023). Finally, we show in Section 4 that for stellar-sized sources to simultaneously meet the composition, energy and density constraints as well as the Hillas criterion, the net closes in on long-duration gamma-ray bursts (L-GRBs).

## 2. Description of the Model

### 2.1. Distribution of Stellar Matter

To model the 3D distribution of galaxies hosting transient UHECR accelerators, we use the near-infrared flux-limited

sample of Biteau (2021) that maps both star formation rate (SFR) and stellar mass ($M_*$) over 90% of the sky. The catalog offers unprecedented cosmography of SFR and $M_*$ within 350 Mpc, including about 400,000 galaxies. This number of entries is an order of magnitude larger than in the Two Micron All-Sky Survey redshift survey (Huchra et al. 2012) used in other recent UHECR studies (e.g., Allard et al. 2022, 2024). Distances are estimated with a 50:50 ratio of spectroscopic and photometric measurements. The catalog also provides corrections for incompleteness with increasing distance and with decreasing Galactic latitude.

Beyond the Local Volume, i.e., at distances greater than 11 Mpc, the incompleteness along the Galactic plane cannot be compensated for. The galaxies above and below the so-called zone of avoidance, which covers the remaining 10% of the sphere, are mirror-symmetrically cloned to fill it, following the approach of Lavaux & Hudson (2011). The cloning procedure has a negligible impact below 50 Mpc, but it allows us to fully reconstruct the contribution of more distant structures, including our supercluster Laniakea whose core region (50–100 Mpc) is close to the zone of avoidance. The catalog of Biteau (2021) thus provides a coherent estimate of the density of matter in the entire 350 Mpc radius volume, with the greatest granularity to date.

The cosmography provided by the catalog is fully consistent with the known distribution of matter in the nearby Universe. Most of the galaxies within a box of 10 Mpc width are in a plane defining the Local Sheet, which is tilted by 8° with respect to the supergalactic plane (McCall 2014). The $M_*$ distribution within a box of 600 Mpc width matches the 3D dark-matter density fields inferred from Cosmicflows-2 (see Hoffman et al. 2017), which models the dark-matter density in a comparable volume using the peculiar velocities of around 8000 galaxies (see Bister & Farrar 2024, and references therein for recent applications to UHECR studies on large angular scales). Finally, the SFR and $M_*$ densities accounting for incompleteness corrections converge beyond 100 Mpc toward the values measured in deep fields (Driver et al. 2018), as illustrated in Figure 1.

### 2.2. Source-emission Model

In our model, each galaxy in the catalog is supposed to harbor UHECR sources, which burst at a given rate. We consider two tracers for the burst rate: SFR and $M_*$. In both cases, the intensity observed on Earth receives contributions from each galaxy in the catalog within 350 Mpc and from a continuous density of galaxies beyond 350 Mpc. By using these two different tracers, we can show the limited influence on our results of the source evolution at high redshift ($z > 0.1$).

The emission rate per unit energy of a single galaxy $i$ is conveniently parameterized as

$$ \dot{Q}_{Ai}(E) = Q_{A}^{0}f(E, Z_{A})S_{i}, \qquad (1) $$

with $S_i$ the ejection rate of UHECRs and $Q_A^0 f(E, Z_A)$ is the number of particles per unit energy. The widely used exponentially suppressed power-law function is adopted for $f(E, Z_A)$, as it provides a summary of the acceleration and ejection processes with only three parameters: the spectral index of protons, the spectral index of heavier nuclei, and the rigidity-dependent maximum energy (see Section 3.1 for details). The emission spectrum is supposed identical at each
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The SFR density (SFRD; left) and stellar-mass density (SMD; right) as a function of distance. A Chabrier initial mass function is adopted for the data and models. The density profiles within a 350 Mpc radius ($z_{\text{cat}}^{\text{max}} \approx 0.08$) are calculated from the catalog of Biteau (2021). The green bands illustrate the 68% confidence interval of deep-field measurements at $0.02 < z < 0.08$ (Driver et al. 2018). The dashed black lines show the cosmic evolution in best agreement with the two tracers (López Fernández et al. 2018) as parameterized and normalized by Koushan et al. (2021). For comparison with the measured profiles, extrapolations of cosmic evolutions at redshifts $z < 0.08$ are represented by gray dotted lines.</figcaption></figure>

<figure><img src="image_3.png" /></figure>

source. Observational data indeed suggest little variation in maximum rigidity from one source to another, to preserve the mixture of elements as the mass composition gets heavier above EeV energies (Ehlert et al. 2023). The ejection rate, on the other hand, is governed by the SFR or $M_*$ of the galaxy

$$S_i = \begin{cases} k \frac{\text{SFR}_i}{c_{N,s}(d_i, b_i)}, \\ \dot{k} \frac{M_{*i}}{c_{N,m}(d_i, b_i)}, \end{cases} \qquad (2)$$

with $k(\dot{k})$ expressed in $M_{\odot}^{-1} (M_{\odot}^{-1} \text{yr}^{-1})$ units and where $c_{N,s}$ ($c_{N,m}$) is the incompleteness correction factor. This correction factor, which depends on the luminosity distance ($d_i$) and Galactic latitude ($b_i$) of the galaxy indexed by $i$, compensates for fainter galaxies below the flux threshold of the near-infrared survey. These faint galaxies are assumed in Equation (2) to follow the same spatial distribution as the brighter galaxies in the catalog. In the following, it will prove useful to express $\dot{Q}_{Ai}(E)$ as

$$\dot{Q}_{Ai}(E) = \frac{p_A(E)}{E} \mathcal{E}_A S_i, \qquad (3)$$

with $\mathcal{E}_A = Q_A^0 \int dE \, \text{Ef}(E, Z_A)$ the total energy emitted per source of UHECRs with mass number A above a minimum rigidity $E/Z_A e$. The minimum rigidity is set here at 0.6 EV by the minimum energy of the observed proton spectrum. This way $p_A(E) = Q_A^0 \text{Ef}(E, Z_A)/\mathcal{E}_A$ is the probability density function for the emitted energy E of an UHECR with mass number A.

Beyond $d_{\text{cat}}^{\text{max}} = 350$ Mpc, that is a redshift $z_{\text{cat}}^{\text{max}} \approx 0.08$, the near-infrared survey used in Biteau (2021) is too shallow to resolve more than half of the matter density in individual galaxies. The ejection rate is then expressed per unit energy and per comoving unit volume as

$$\dot{q}_A(E, z) = q_A^0 f(E, Z_A) \begin{cases} k \dot{\rho}_{\text{SFR}}(z). \\ \dot{k} \rho_{M_*}(z). \end{cases} \qquad (4)$$

Here, $\dot{\rho}_{\text{SFR}}(z)$ and $\rho_{M_*}(z)$ are the SFR density (SFRD) and stellar-mass density (SMD), respectively. We adopt the cosmic

evolution of the SFRD, $\dot{\rho}_{\text{SFR}}(z)$, as provided by López Fernández et al. (2018) and normalized by Koushan et al. (2021). Of all the models tested by Koushan et al. (2021), the parameterization adopted here provides the best agreement with both SFR and $M_*$ measurements at redshifts $z \le z_{\text{cat}}^{\text{max}}$ (Driver et al. 2018; Biteau 2021). The model for the cosmic evolution of the SMD, $\rho_{M_*}(z)$, is obtained by integrating the SFRD over lookback time, the role of which is played by redshift

$$\rho_{M_*}(z) = (1-R) \int_z^\infty dz' \dot{\rho}_{\text{SFR}}(z') \left| \frac{dt}{dz'} \right|, \qquad (5)$$

with $R=0.41$ the return fraction for a Chabrier initial mass function (Madau & Dickinson 2014).$^4$ This fraction corresponds to the proportion of the mass of each generation of stars that is returned to the interstellar and intergalactic media. The term $|dt/dz|$ follows from the Lambda cold dark matter (ΛCDM) cosmological parameters

$$\left| \frac{dt}{dz} \right| = \frac{1}{H_0 (1+z) \sqrt{(1+z)^3 \Omega_m^0 + \Omega_\Lambda^0}}, \qquad (6)$$

with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m^0 = 0.3$, and $\Omega_\Lambda^0 = 0.7$.

## 2.3. Intensities on Earth

The all-particle energy spectrum observed on Earth, $J(E)$ in units of $\text{eV}^{-1} \text{km}^{-2} \text{yr}^{-1} \text{sr}^{-1}$, receives contributions from the foreground galaxies in the catalog and from an isotropic background at $z > z_{\text{cat}}^{\text{max}}$

$$J(E) = \sum_A \left[ J_A^{\text{bkg}}(E) + \frac{1}{4\pi} \sum_i \int dn \, J_A^i(E, n) \right], \qquad (7)$$

where the outer sum is over each observed species of mass number A and the inner sum is over the galaxies in the catalog. We explain below the scenario tracing the UHECR sources by the SFR of galaxies. The scenario based on stellar mass is obtained by substituting $k\rho_{\text{SFR}}(z) (k\text{SFR}/c_{N,s})$ with $\dot{k}\rho_{M_*}(z) (\dot{k}M_* / c_{N,m})$.

<footnote>$^4$ The choice of initial mass function does not affect our results on the absolute density of the UHECR burst rate, expressed in units of $\text{Gpc}^{-3} \text{yr}^{-1}$.</footnote>
Up to $z = z_{\text{cat}}^{\text{max}}$, each single galaxy contributes

$$ J_A^i(E, \mathbf{n}) = \frac{F(\mathbf{n}, \mathbf{n}_i; \Theta)}{4\pi d_i^2} \sum_{A'} \int dE' \dot{Q}_{A'}(E') \frac{d\eta_{AA'}(E, E', z)}{dE} \\ = \frac{SFR_i}{4\pi d_i^2 c_{N,s}(d_i, b_i)} F(\mathbf{n}, \mathbf{n}_i; \Theta) \sum_{A'} k\mathcal{E}_{A'} \\ \times \int \frac{dE'}{E'} p_{A'}(E') \frac{d\eta_{AA'}(E, E', z)}{dE}, \qquad (8) $$

to the observed spectrum for a given species. The factor $\mathcal{E}_{A'}$ is defined here as the energy radiated isotropically in each emitted nuclear species $A'$. The energy losses and spallation processes are described by $\eta_{AA'}(E, E', z)$, which is the fraction of particles detected on Earth with energy $E$ and mass number $A$ from parent particles emitted by the sources with energies $E' \gg E$ and mass numbers $A' \gg A$. In practice, $\eta_{AA'}(E, E', z)$ is calculated by Monte Carlo using the SimProp package (Aloisio et al. 2017).<fnref n="5" /> Finally, $F(\mathbf{n}, \mathbf{n}_i; \Theta)$ represents the probability density function of the arrival direction $\mathbf{n}$ given the direction of the galaxy $\mathbf{n}_i$. We use a von Mises-Fisher distribution (equivalent to a Gaussian distribution on a sphere) centered on $\mathbf{n}_i$ with an adaptable parameter $\Theta$, which is aimed at modeling in an effective way the magnetic deflections experienced by UHECRs throughout their propagation. This choice for $F$ is based on the assumption that some directionality of the sources is preserved in the observed sky maps due to the high rigidity of the particles, a reasonable assumption based on what is known on the strength of the magnetic fields experienced. A comprehensive treatment of the deflections would require precise knowledge of the Galactic magnetic field and of the filamentary structure and strength of the Local Sheet magnetic field, both of which are expected to imprint random-type deflections from magnetic turbulence as well as mass-spectrometric-type deflections from large-scale fields. Despite recent progresses, such knowledge is still elusive. Our above choice allows us to emulate deflections to first order with a single model parameter, $\Theta$. Some comparisons to simulations with models of magnetic fields are further discussed in the Appendix.

For $z > z_{\text{cat}}^{\text{max}}$, the energy spectrum $J_A^{\text{bkg}}(E)$ observed at present time for a given species $A$ is modeled from the background density of sources as

$$ J_A^{\text{bkg}}(E) = \frac{c}{4\pi} \int dz \left| \frac{dt}{dz} \right| \sum_{A'} \int dE' \dot{q}_{A'}(E', z) \frac{d\eta_{AA'}(E, E', z)}{dE} \\ = \frac{c}{4\pi} \sum_{A'} k \mathcal{E}_{A'} \int dz \left| \frac{dt}{dz} \right| \dot{p}_{\text{SFR}}(z) \\ \times \int \frac{dE'}{E'} p_{A'}(E') \frac{d\eta_{AA'}(E, E', z)}{dE}. \qquad (9) $$

Here, the factor $c/4\pi$ converts the local density of UHECRs per unit energy into a flux per steradian (intensity) and per unit energy.

This background contribution is assumed isotropic, a reasonable assumption beyond a few hundred megaparsecs (Biteau 2021).

## 2.4. Time Delays and Magnetic Deflections

In the transient scenario, sources are visible by an observer only during a finite duration $\Delta\tau$. Due to the magnetic wandering of particles en route to the Earth, the duration window depends on the magnetic rigidities of the particles, the distance from the source, and the magnetic-field properties. For an observation duration much smaller than $\Delta\tau$, the mean number of relic bursts that contribute to the energy spectrum observed today on Earth reads as

$$ \lambda_i = S_i \Delta\tau. \qquad (10) $$

The stochastic nature of any transient scenario can thus be captured in Equation (8) by randomly drawing the number of bursts in each galaxy according to a Poisson distribution of parameter $\lambda_i$. The median scenario, also known as the Asimov model, can serve as a benchmark for data comparison. The median $\mu_{1/2}$ of the Poisson draws is conveniently approximated here by

$$ \mu_{1/2} = \begin{cases} \lambda & \text{for } \lambda > \lambda_{\text{th}}, \\ 0 & \text{otherwise,} \end{cases} \qquad (11) $$

with $\lambda_{\text{th}} \approx \ln 2$. We have verified that this approximation for the Poisson process provides the expected median over 100 random draws of sky maps, while also considerably simplifying the computation.<fnref n="6" />

The image on Earth of a burst of duration $\delta t$ is visible during a period $\Delta\tau \gg \delta t$, which can be shown to scale with the amplitude of the magnetic field, $B$, the particle rigidity, $R$, the source distance, $d$, and the coherence length of the field, $\lambda_B$, as (e.g., Achterberg et al. 1998)

$$ \begin{aligned} \Delta\tau &= \frac{\sqrt{2} \lambda_B}{9c} \left(\frac{cBd}{R}\right)^2 \\ &= 70 \text{ kyr} \times \left(\frac{B}{10 \text{ nG}}\right)^2 \left(\frac{R}{5 \text{ EV}}\right)^{-2} \left(\frac{d}{2 \text{ Mpc}}\right)^2 \left(\frac{\lambda_B}{10 \text{ kpc}}\right). \end{aligned} \qquad (12) $$

Depending on the intervening magnetic fields, the temporal spread $\Delta\tau$ can be substantial and make a transient event quasi persistent for an observer far enough away from the source (see, e.g., Sigl 2001 for a review). The magnetic fields of interest that may enter into Equation (12) can be found in four environments: the surroundings of the source (see Section 2.5), the extragalactic medium, the Local Sheet, and the Milky Way.

Upper limits on extragalactic magnetic fields in cosmic voids are conservatively set to the nG level from rotation measurements (Pshirkov et al. 2016; O'Sullivan et al. 2020). Yet, they can be refined down to 10–50 pG for magnetic fields of primordial origin that would affect CMB anisotropies (Jedamzik & Saveliev 2019). On the flip side, lower limits at the fG level have also been derived from the nonobservation in the GeV range of gamma-ray cascades from TeV blazars (Neronov & Vovk 2010; Tavecchio et al. 2010; Ackermann et al. 2018;

<footnote n="5">⁵ We use SimProp v2r4 with the model from Gilmore et al. (2012) for the extragalactic background light and the photo-disintegration cross sections from Koning et al. (2005) and Koning & Rochman (2012), setup with $M = 4$ (see Aloisio et al. 2014).</footnote>

<footnote n="6">⁶ The careful reader will note that this approximation violates the additive property of Poisson variables. Verification using random draws of sky maps is a necessary step in the study.</footnote>
Aharonian et al. 2023). Therefore, we adopt here a conservative value of 10 pG for a coherence length of 1 Mpc, as suggested by CMB constraints for a primordial origin of extragalactic magnetic fields and by cosmological magnetohydrodynamics simulations for an astrophysical origin (Vazza et al. 2017). With such characteristics, extragalactic magnetic fields in voids induce a time spread of ≈70,000 yr for UHECRs at a magnetic rigidity of 5 EV from sources at 200 Mpc.

The Local Sheet refers to the planar structure of the cosmic web that includes the Local Group of galaxies and the ring of surrounding galaxies known as the Council of Giants. According to McCall (2014), the Local Sheet can be described as a disk of diameter ~10 Mpc and thickness ~0.5 Mpc. About 1 Mpc from the center of the sheet are Andromeda and the Milky Way. These two galaxies are at the heart of the Local Group, which can be described as a disk nearly coplanar with the Local Sheet, with a diameter of ~5 Mpc and a thickness of ~1.5 Mpc. The plane of the Milky Way is nearly perpendicular to that of the Local Sheet, as might be expected if the spin of galaxies is determined by the tidal field during the formation of the sheet through gravitational collapse (see, e.g., Aragon-Calvo et al. 2023, and references therein).

The magnetic field permeating the Local Sheet is largely underconstrained, as is the distribution of magnetic fields in the cosmic web. The field at the nodes of the web is constrained in the brightest galaxy clusters. The warm-hot plasma that resides in the intracluster medium embeds a turbulent magnetic field with a coherence length of ≈10 kpc and an amplitude of the order of a few µG, whose fluctuations are well represented by a Kolgomorov power spectrum (e.g., Bonafede et al. 2010, for the Coma cluster). The field in the structures connecting the nodes has been detected by radio and X-ray stacking of galaxy filaments (Vernstrom et al. 2021). The authors estimate a magnetic-field amplitude in filaments of the order of 30 nG, a value that is confirmed by Faraday-rotation measurements (Carretti et al. 2022). The latter analysis suggests a predominantly structured field in the filaments, with a turbulent component of amplitude estimated at 3–6 nG. The field in the sheets between the filaments is unconstrained, both in observations and in cosmological simulations. We assume a warm-hot plasma with properties similar to those of clusters, i.e., a coherence length $\lambda_B \approx 10$ kpc. The turbulent field is presumed to have an amplitude comparable to that observed in filaments. Given the paucity of constraints on cosmic magnetism, we consider a priori that plausible values are between 1 and 100 nG. We adopt these values to model the Local Sheet field as Kolmogorov turbulence extended over the radius of the Local Group, 1–3 Mpc. Such radii correspond to the underdensity range in Figure 1 (see also McCall 2014, for a detailed discussion of the size and geometry of the Local Group). The geometric description of the magnetic field of the Local Sheet as a central sphere encompassing the Local Group is the simplest approximation that can be tested at this stage, since the properties of the plasma at the heart of the Local Sheet are still poorly understood both theoretically and observationally.

With a radius of 2 Mpc, a magnetic field of 10 nG strength, and a 10 kpc coherence length injected into Equation (12), the time spread of UHECRs at a magnetic rigidity of 5 EV accumulated in the Local Sheet amounts to ≈70,000 yr for sources beyond 2 Mpc. Consequently, the Local Sheet field is expected to show a larger contribution to $\Delta\tau$ than the extragalactic magnetic field in voids for sources within the

<table><caption>**Table 1**
Default Parameters for the Turbulent Magnetic Fields Considered in the Model</caption><thead><tr><th></th><th>r<sub>max</sub><br>(Mpc)</th><th>B<br>(nG)</th><th>&lambda;<sub>B</sub><br>(kpc)</th></tr></thead><tbody><tr><td>Milky Way</td><td>10<sup>-2</sup></td><td>10<sup>3</sup></td><td>10<sup>-1</sup></td></tr><tr><td>Local Group/Sheet</td><td>2</td><td>10</td><td>10</td></tr><tr><td>Voids</td><td>&infin;</td><td>10<sup>-2</sup></td><td>10<sup>3</sup></td></tr></tbody></table>

**Note.** For each field, uniform turbulence is assumed to fill a spherical region of radius $r_{max}$ (column (2)), with an amplitude $B$ (column (3)) and coherence length $\lambda_B$ (column (4)).

GZK horizon. This structure also contributes to deflect particles following (see, e.g., Achterberg et al. 1998)

$$ \begin{aligned} \Delta\theta &= \frac{2 c B \sqrt{\lambda_B d}}{3 R} \\ &= 10^{\circ} \times \left(\frac{B}{10 \text{ nG}}\right) \left(\frac{R}{5 \text{ EV}}\right)^{-1} \left(\frac{d}{2 \text{ Mpc}}\right)^{1/2} \left(\frac{\lambda_B}{10 \text{ kpc}}\right)^{1/2} \end{aligned} \quad (13) $$

The spatial extension of the Galactic magnetic field is too limited to imprint sizeable time spreads for extragalactic sources. However, with a strength of several microgauss in the disk, as estimated from polarization and Faraday-rotation measurements, the Galactic magnetic field is sufficiently intense to contribute to UHECR angular deflections. We account for the corresponding contributions in the Appendix with the model of Jansson & Farrar (2012). To explore the parameter space more feasibly, we approximate the Galactic magnetic field by Kolmogorov turbulence extending to 10 kpc, with an amplitude of 1 µG and a coherence length of 100 pc. For these values, an UHECR emitted at a distance $d > 10$ kpc with a rigidity of 5 EV undergoes in the Milky Way a time delay of ≈200 yr and an angular deviation of ≈7°, which add quadratically to those induced by the Local Sheet. A summary of the default parameters used for the turbulent components of the magnetic fields relevant to this study is provided in Table 1.

The effect of the delay in Equation (12) is conveniently implemented by noting that Equations (10) and (11) impose a maximum rigidity for each galaxy of ejection rate $S_i$ and distance $d_i$.<fnref n="7" /> Above this rigidity, UHECR bursts are too short and infrequent to be visible in more than half of all realizations. The median sky maps simulated above an energy threshold take into account the angular spread using a von Mises-Fisher distribution, as described in Equation (8). The angular-spread parameter is determined using Equation (13) for $R = \bar{R}$, where $\bar{R}$ is the average rigidity of the UHECRs above the energy threshold, i.e., $\Theta = \Delta\theta(\bar{R})$.

## 2.5. Galaxy-cluster Extinction

The impact of the source environment on the observed angular extent of UHECR excesses can be neglected for a distance from the source significantly greater than the size of

<footnote n="7">⁷ The conservation of the Lorentz factor in the photodissociation process allows the observed rigidity to be approximated by the emitted rigidity of UHECR nuclei. This approximation is limited here to the implementation of time delays, including the effect of the clusters in Section 2.5. All energy losses are otherwise fully taken into account as shown in Equations (8) and (9).</footnote>
the environment. The time delay caused by a source environment less magnetized than the Local Sheet is effectively taken into account through the threshold of the Poisson process in Equation (11). Only an environment with a turbulent magnetic field well above the range adopted for the Local Sheet (1–100 nG) could affect the simulated sky maps. This is precisely the case for galaxy clusters, which confine UHECRs long enough to act as effective calorimeters (Dolag et al. 2009; Kotera et al. 2009; Harari et al. 2016; Fang & Murase 2018). We model here the effect of galaxy clusters following Condorelli et al. (2023), whose approach is summarized below.

In Condorelli et al. (2023), the hot, tenuous gas of the intracluster medium is described under the self-similar assumption for each galaxy cluster, so that its radial pressure profile can be described using only its mass and redshift (Kaiser 1986). The magnetic field is scaled to the thermal components of the intracluster medium, assuming that the magnetic energy density is proportional to the thermal energy density, with a relation calibrated using the Perseus cluster ((l, b) = (150°6, -13°3), d ≃ 70 Mpc, M ≃ 5.8 × 10¹⁴ M☉, and R₅₀₀ ≃ 1.3 Mpc, from Urban et al. 2014). With these ingredients in hand, the propagation of UHECRs along the gas-density and magnetic-field profiles can be simulated inside any galaxy cluster. A modified version of SimProp, which includes hadronic processes and magnetic-field confinement, has enabled Condorelli et al. (2023) to construct generic transparency functions of clusters as a function of UHECR rigidity for protons on the one hand, and nuclei on the other.

The brightest X-ray galaxy cluster in the GZK horizon is the Virgo cluster ((l, b) = (279°7, 74°5), d ≃ 15 Mpc, M ≃ 1.2 × 10¹⁴ M☉, and R₅₀₀ ≃ 0.8 Mpc, from Planck Collaboration et al. 2016), followed by the Perseus cluster, which is 1.5–2 times fainter in X-rays. These clusters are opaque to both nuclei and protons in the rigidity range of interest. We apply the transparency function of Condorelli et al. (2023) to all galaxies behind disks of radius 3 R₅₀₀ centered on Virgo and Perseus. The next cluster is Coma, which is ≃ 5 times fainter than Virgo in X-rays. We have verified that implementing the UHECR opacity of Coma or fainter clusters has no impact on the reconstructed sky maps in the energy range of interest. We do note, however, that mass-composition anisotropies might be expected in the energy range where protons and heavier nuclei mix (near the ankle), due to the different filtering that these species undergo in the traversed intracluster media. We leave the exploration of the UHECR sky around the ankle to future studies that may elaborate on this framework.

# 3. Constraints on Burst-rate Density from Ultrahigh-energy Cosmic-Ray Observations

## 3.1. Best Fit of the Spectrum and Mass-composition Data

The spectrum and mass-composition data used in this work are those obtained at the Pierre Auger Observatory (Aab et al. 2015), which provides the largest exposure to date to UHECRs. The Pierre Auger Observatory is a hybrid system that detects the extensive air showers induced by the collisions of UHECRs with atmospheric molecules. We use the all-particle energy spectrum inferred from these data (Abreu et al. 2021) and a proxy of the primary mass of the particles (Aab et al. 2014b; Fitoussi et al. 2023), i.e., the slant depth of maximum of shower development, Xmax. Using

hadronic-interaction generators to model the development of the showers, the Xmax distributions allow the energy-dependent mass composition to be deduced on a statistical basis. As in Tkachenko et al. (2023), we used the two up-to-date hadronic-interaction generators that best describe the Xmax distributions, namely EPOS-LHC (Pierog et al. 2015) and Sibyll2.3d (Riehn et al. 2020). We adopt the latter model in our reference configuration.

The best fit of the spectrum and mass-composition data is obtained following Luce et al. (2022), with a source evolution following either the SFRD or the SMD shown in Figure 1. With more details than given in Section 2.2, we generically model the ejection rate per comoving unit volume and per unit energy of nucleons as

$$ \dot{Q}_{\mathrm{p}}(E) = \dot{Q}_{\mathrm{p}}^{0} \left( \frac{E}{Z_{\mathrm{p}}E_{0}} \right)^{-\gamma_{\mathrm{p}}} f_{\mathrm{supp}}(E, Z_{\mathrm{p}}), \quad (14) $$

with Zp = 1, and a single reference ejection rate Qp0, spectral index γp, and suppression function fsupp(E, Zp) for both escaping protons and protons from neutron decay. Here, E0 is arbitrarily set to 1 eV. The ejection rate of nuclei with mass number Ai is also generically modeled as

$$ \dot{Q}_{A_i}(E) = \dot{Q}_{A_i}^{0} \left( \frac{E}{Z_{A_i}E_0} \right)^{-\gamma_{A_i}} f_{\text{supp}}(E, Z_{A_i}), \quad (15) $$

which describes a single spectral index γA and four independent reference ejection rates QA0 for helium, nitrogen, silicon, and iron. The suppression function adopted both for nucleons and nuclei is taken as

$$ f_{\text{supp}}(E, Z_A) = \begin{cases} 1 & \text{if } E \le E_{\text{max}}^{Z_A}, \\ \exp(1 - E/E_{\text{max}}^{Z_A}) & \text{otherwise.} \end{cases} \quad (16) $$

The maximum acceleration energy is assumed to be proportional to the electric charge of each nucleus, Emax = ZA Emax, with a single free parameter Emax shared by all species. To connect properties at the sources with observables at Earth, the extragalactic propagation is carried out with SimProp. All the events generated are stored in a normalized 5D array T(Edet, Adet | z, Einj, Ainj). In this way, the tensor encodes the probability that a particle with energy Einj and mass Ainj emitted at a redshift z is detected with an energy Edet and a mass Adet. After taking into account the possible interactions of UHECRs with photons of the CMB and of the extragalactic background light during extragalactic propagation, the flux at Earth can be calculated and then compared to the experimental data. This comparison is made using a deviance, D, similar to that presented in Equation (6) of Luce et al. (2022). The first deviance term is defined here as a χ² including the points of the proton spectrum at energies below the ankle, the all-particle spectrum above the ankle, and the mean and standard deviation of Xmax above the ankle (red and black markers, respectively, in Figure 2). The second deviance term accounts for the systematic uncertainty on the mean Xmax values, σsys. = 8 g cm⁻², so that D = χ² + (δXmodel/σsys.)², where δXmodel is an energy-independent shift applied to the model.

We present the results obtained with a source evolution following the SFRD in Figure 2. Consistent to that found in
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The energy spectrum (left) and $X_{max}$ moments (right) of UHECRs as a function of energy. The data points are from the latest data sets of the Pierre Auger Observatory (Abreu et al. 2021; Fitoussi et al. 2023), with a proton fraction based on the Sibyll2.3d hadronic-interaction model (see Tkachenko et al. 2023). The best-fit model, including a shift in $X_{max}$ scale by $\delta X_{model} = (1.20 \pm 0.18) \times \sigma_{sys.} = 9.6 \pm 1.4$ g cm$^{-2}$, is shown as a brown solid line in its range of applicability, and as a brown dashed line at lower energies.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

other studies, the value of $E_{max}$ is determined primarily by the drop in the proton component at $E_{max} = (1.9 \pm 0.1)$ EeV. The spectral index of the nuclei, $\gamma_A = -0.36 \pm 0.21$, is in turn determined by the increase of the average mass with energy, which is almost mono-elemental. The best-fit solution therefore imposes a hard index for nuclei so that the contribution of each element mixes as little as possible. Protons are present in an energy range where a mixture of elements is required to model the all-particle spectrum and the composition. The best-fit value of $\gamma_p = 2.6 \pm 0.7$ is much softer than that of $\gamma_A$. The spectral index of the proton component is affected by a large uncertainty. Indeed, in a scenario where the density of sources is larger in the past, a nonnegligible fraction of the protons observed today originate from the interaction of heavier nuclei during their extragalactic propagation. Thus, for a total emissivity per unit SFR of $(4.8 \pm 0.2) \times 10^{46}$ erg $M_{\odot}^{-1}$, only $(19 \pm 4)$% of the energy is emitted as protons, while the proportions of energy emitted as He, CNO, Si, and Fe are $(14 \pm 3)$%, $(48 \pm 4)$%, $(15 \pm 3)$%, and $(4 \pm 1)$%, respectively.

Note that the deviance of the fit, $D = 54.3$ for 36 degrees of freedom, is significantly improved by shifting the $X_{max}$ scale of the model for UHECR protons nuclei by $\delta X_{model} = (1.20 \pm 0.18) \times \sigma_{sys.}$ This would correspond to a data shift by $-(1.20 \pm 0.18) \times \sigma_{sys.}$ which is in line with values ranging from $-0.9 \times \sigma_{sys.}$ to $-0.5 \times \sigma_{sys.}$ as identified already in Aab et al. (2017a). Systematic uncertainties in the $X_{max}$ scale, $\delta X_{model}$, should also impact the proton spectrum below the ankle, but this effect is difficult to take into account with the public data, which are provided for a zero value of this nuisance parameter only. The systematic uncertainty of $\sigma_E/E \simeq 14\%$ in the energy scale has a significant impact on the all-particle and proton spectra. A shift of the model by $+1\sigma_E$ leads to a lower emissivity $(2.4 \pm 0.2) \times 10^{46}$ erg $M_{\odot}^{-1}$ and a better deviance ($D = 46.0$). A shift by $-1\sigma_E$ leads to higher emissivity $(9.7 \pm 0.4) \times 10^{46}$ erg $M_{\odot}^{-1}$ and a worse deviance ($D = 61.4$).

For a source evolution following the SMD, the emissivity per unit $M_*$ is estimated to $(4.2 \pm 0.2) \times 10^{36}$ erg $M_{\odot}^{-1}$ yr$^{-1}$, with a cutoff energy at $2.6 \pm 0.1$ EeV, a proton index of $3.0 \pm 0.2$, and a nuclear index of $0.6 \pm 0.1$. The emitted proton component, which saturates the observed fraction, carries $(54 \pm 4)$% of the energy, $(32 \pm 4)$% being carried by CNO and less than 5% carried by each of the other species. The deviance

in the SMD scenario reaches 68.3 units for a shift in the $X_{max}$ scale of the model by $(1.5 \pm 0.2) \times \sigma_{sys.}$ These fits of the model to the data can be easily reproduced using the public version of the analysis code, which we make available in Zenodo (see Condorelli et al. 2024).

## 3.2. Full-sky Observations above $\simeq 40$ EeV

The Pierre Auger Observatory, covering $3000$ km$^2$ in the province of Mendoza (Argentina), has since 2004 accumulated a total exposure of $\simeq 135,000$ km$^2$ sr yr over nearly the entire UHECR sky, except in the northernmost quarter. In parallel, the Telescope Array (Abu-Zayyad et al. 2013), covering $700$ km$^2$ in Utah (USA), has provided since 2008 a view of the northern UHECR sky with a total exposure of $\simeq 17,500$ km$^2$ sr yr. A series of intercollaborative efforts has resulted in the full-sky map of arrival directions shown in Figure 3 in Galactic coordinates (Abdul Halim et al. 2023), with top-hat smoothing on an angular scale $\Psi = 25^\circ$ (corresponding to a Gaussian angular scale $\Theta_{obs} \approx \Psi/1.59 \simeq 16^\circ$; see the discussion in Abreu et al. 2022). Although the techniques for assigning energies to events are nearly the same at both observatories, there are differences in the way primary energies are derived. To produce the sky map, the energy estimators are cross-calibrated by ensuring that the intensities observed at each observatory are identical in the area common to both fields of view. On account of this, the sky map is built with events at energies above 38 EeV (48 EeV) in terms of the energy scale of the Pierre Auger Observatory (Telescope Array). For convenience, we choose to refer to $\simeq 40$ EeV for the energy threshold hereafter. Note that, as accumulated exposure in the Northern Hemisphere is around 8 times lower than in the Southern Hemisphere, the fluxes reached by the excesses are prone to greater statistical uncertainties on the northernmost quarter of the sky, i.e., on the left quarter of the map in Galactic coordinates in Figure 3. These uncertainties are illustrated by the Poisson fluctuations visible at an angular scale of $\simeq 5^\circ$, which corresponds to the typical spacing between each of the $\sim 2000$ UHECRs detected over the entire sphere at $E \gtrsim 40$ EeV.

Due to the GZK effect, the isotropic background of UHECR arrival directions caused by sources distributed throughout the distant Universe is suppressed. The overdensities hinted at in Figure 3 could therefore point to dominating foreground sources. However, a blind search for excesses does not reveal
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Flux map of UHECRs at energies above 38 EeV (48 EeV) as reconstructed from the Pierre Auger Observatory (Telescope Array). The sky map is smoothed using a top-hat function of radius $\Psi = 25^{\circ}$. The gray line shows the supergalactic plane. Adapted from Abdul Halim et al. (2023).</figcaption></figure>

any statistically significant indications for anisotropy (di Matteo et al. 2020; Abreu et al. 2022). Indeed, by not specifying a priori the targeted regions of the sky where the excesses are searched for, nor the angular-window radius or the energy threshold, the analysis sensitivity suffers from the numerous trials performed.

Yet, even without compelling proof for discrete sources, a correlation between UHECR arrival directions and the flux pattern of a class of astrophysical objects might reveal earlier an anisotropy that would trace the sources. The best evidence is currently provided by a catalog of bright star-forming or starburst galaxies (Aab et al. 2018; Abreu et al. 2022) with a confidence level of $4.6\sigma$ when considering full-sky maps (Abdul Halim et al. 2023). The analysis assesses the clustering of UHECR arrival directions around the predefined astro-physical source positions, in proportion to the nonthermal photon luminosities of the host galaxies. The similarities between the model and observations, best captured on an angular scale of $\Theta_{\text{obs}} = (15^{+5}_{-3})^{\circ}$, are due to overdensities in the direction of Local Sheet galaxies. However, the low signal strength, $(12^{+5}_{-3})\%$, precludes the conclusion that UHECRs come exclusively from these bright nearby galaxies.

Interestingly, the brightest star-forming galaxies known as starburst galaxies turn out to be responsible for ≈15% of the total SFR for redshifts $z < 2$ (Rodighiero et al. 2011; Sargent et al. 2012). If UHECR sources are traced by the SFR (or $M_*$), then the signal intensity could suggest that the sources are hosted by every star-forming galaxy in the Universe. With their higher SFR per unit $M_*$, starburst galaxies would merely have a higher probability of hosting UHECR transient sources than main-sequence galaxies with similar $M_*$. The sky maps constructed below allow us to explore the rate of bursts in each galaxy relative to its SFR (or $M_*$) that could, given the observed angular scale of the correlation, accommodate the observed signal strength.

### 3.3. Constraints on k from Arrival Directions

As shown in Equation (10), the number of bursts in each galaxy can be drawn at random according to a Poisson process of parameter $\lambda_i$. For a set of magnetic fields, the Poisson parameter is governed by the rate of bursts $k$ relative to the SFR (or $M_*$) of the host galaxy (see Equation (2)). For a set of parameters fitted to the energy spectrum and mass-composition data (see Section 3.1), the intensity of UHECRs observed on Earth in each direction above ≈40 EeV is then estimated from Equations (8) and (9). The results are illustrated, for a given $k$ value, by means of the median intensity in each direction. These median sky maps better reflect the stochastic nature of the transient scenario than the mean sky maps, which would just depict the results obtained by weighting each galaxy by $\lambda_i$.

The value of $k$ affects the sky maps as follows. For $k$ values such that the emission rates are large ($\lambda_i > 1$), each galaxy emits almost persistently. The resulting burst rate is large enough to imply dominant contributions from galaxies of the Local Group, in particular Andromeda ($d \approx 750$ kpc), the Magellanic clouds ($d \approx 50–60$ kpc), and the Milky Way. For lower values of $k$, galaxies in the Local Group are filtered out. This is because the time spread induced by the Local Sheet magnetic field is not sufficient to make the contribution of nearby galaxies stationary: their burst rate is so low that the observer is likely to be between two bursts. For even lower values of $k$, the dimmest galaxies in the Local Sheet and beyond disappear due to their low burst rate.

The transient emission from our Galaxy requires special treatment due to its proximity. The sky maps can be evaluated by concentrating all Galactic matter in its bulge, centered on Sgr A*, using the SFR and $M_*$ from Licquia & Newman (2015) normalized to the Chabrier initial mass function. Alternatively, the sky maps can be evaluated omitting the matter contained in the Milky Way, which yields more pessimistic constraints. These assumptions encompass the results expected for a three-dimensional distribution of Galactic matter, the treatment of
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Contrast between maximum brightness in the Cen region and the Southern Hemisphere average. The 68% confidence region from the Pierre Auger Observatory is shown with gray lines (Abreu et al. 2022). The parameter space is explored with a scan of the burst rate, *k* (x-axis), and of the magnetic field in the Local Sheet, *B*<sub>LS</sub> (color bar). The solid red line delimits the parameter space for which localized excesses are observed in the Cen and UMa regions, for a conservative approach excluding the Milky Way from the model. The inclusion of transient emission from our Galaxy in the model, as a point source colocalized with Sgr A*, leads to the allowed region delimited by the gold dotted line.</figcaption></figure>

which within a transient model goes beyond the scope of
this work.

We explore the range of *k* values compatible with the data using the two most prominent excesses identified by the Pierre Auger and Telescope Array collaborations, i.e., the excesses in the Centaurus (Cen) and Ursa Major (UMa) regions, around the directions (*l*, *b*) = (305°, 16°) and (*l*, *b*) = (176°, 45°) in Galactic coordinates. These two excesses are responsible for the evidence of anisotropy above ≈40 EeV (Abdul Halim et al. 2023). The implications of our model for a potential excess in the Andromeda region (eponymous galaxy located at (*l*, *b*) = (121°2, -21°6)) will be discussed.

For a given value of the Local Sheet turbulent magnetic field and of the burst rate, we build the median sky map and search for the local maximum of the flux within a radius of 40° around the Cen excess and that around the UMa excess, as localized by the Pierre Auger and Telescope Array collaborations, respectively. The relatively large value of the search radius (40°) effectively takes into account the displacement of maxima expected as a result of coherent magnetic fields, but whose exact direction remains dependent on the model (Erdmann et al. 2016). If any of the local maxima fall at the edge of its search region, the model is considered to lack one of the hot spots hinted at by the data, and the parameter set is rejected.

The parameter space that allows for UHECR emission traced by the SFR is shown in Figure 4. In an SFR ($M_*$) model, the magnetic field of the Local Sheet is between 1 and 20 nG (0.5–10 nG) for a burst rate $k \in [8 \times 10^{-6}; 2 \times 10^{-3}] M_{\odot}^{-1}$ ($\dot{k} \in [6 \times 10^{-16}; 4 \times 10^{-13}] M_{\odot}^{-1}$ yr$^{-1}$). Including the Milky Way in the model brings the upper bound of the allowed rate range down to $k \le 3 \times 10^{-4} M_{\odot}^{-1}$ ($\dot{k} \le 2 \times 10^{-14} M_{\odot}^{-1}$ yr$^{-1}$). Added to the effect of the Galactic magnetic field, the deflections induced by a 10 nG field in the Local Sheet lead to an angular spread $\Theta \approx 18^\circ$, in agreement with observations ($\Theta_{\text{obs}} = (15^{+5}_{-3})^\circ$; Abdul Halim et al. 2023).

In the suppression region of the UHECR spectrum, given the steep decrease of the number of events with energy, the flux value derived for each hot spot is highly dependent on the energy threshold considered. Rather than the flux, we provide instead the contrast for the hot spot measured with the best accuracy. Data from the Pierre Auger Observatory yield a Cen hot spot flux reaching $\Phi_{\text{Cen}}(>40 \text{ EeV}) = (15.9 \pm 1.3) \times 10^{-3} \text{ km}^{-2} \text{ yr}^{-1} \text{ sr}^{-1}$ in a top-hat region of radius $\Psi = 25^\circ$ (Abreu et al. 2022). This maximum flux reconstructed by the Pierre Auger Observatory can be compared with the average value obtained in the Southern Hemisphere, $\Phi_{\text{SH}}(>40 \text{ EeV}) = (11.3 \pm 0.4) \times 10^{-3} \text{ km}^{-2} \text{ yr}^{-1} \text{ sr}^{-1}$, at decl. $\delta < 45^\circ$. The 68% confidence region on the Cen contrast, $C_{\text{Cen}} = \Phi_{\text{Cen}}/\Phi_{\text{SH}} = 1.41 \pm 0.13$, is shown with gray lines in Figure 4. As illustrated by the colored region in Figure 4, the SFR scenario predicts contrast values of 1.3–2 that are compatible with the observations in the Southern Hemisphere. A somewhat narrower range of contrast values, 1.5–1.85, is obtained in the $M_*$ scenario, which is also compatible with the observations. In both scenarios, comparable contrasts are predicted for the UMa and Cen hot spots. Despite the good agreement between the expected and observed contrasts, we do not use this observable to constrain the parameter space, either for the magnetic field in the Local Sheet or for the UHECR burst rate. The flux in these regions is indeed subject to the effects of deflection and (de) magnification by the coherent components of the traversed magnetic fields (see the Appendix).

The contrast obtained for the third hot spot in the Andromeda region (see Figure 5, where the location of the eponymous galaxy is marked) depends strongly on the region of parameter space considered. This statement is illustrated by Figure 5, where the sky model is explored for different values of the burst rate and for a fixed amplitude of the turbulent magnetic field in the Local Sheet. For the highest burst rates ($k \ge 2 \times 10^{-4} M_{\odot}^{-1}$), the Andromeda and Triangulum galaxies produce an excess brighter than those from the Local Sheet galaxies and from our supercluster Laniakea. For $k = 1 \times 10^{-4} M_{\odot}^{-1}$, illustrated by the model map in Figure 5, which appears to match the observed map in Figure 3, the catalog up to $d_{\text{cat}}^{\text{max}} = 350 \text{ Mpc}$ accounts for ≈97% of the flux above 40 EeV, with only ≈3% of the flux coming from the sphere of radius 1 Mpc that includes Andromeda and Triangulum, ≈20% from the 1 to 10 Mpc shell that includes the galaxies of the Council of Giants, and ≈50% from the shell between 10 and 100 Mpc that contains our supercluster Laniakea. Laniakea contributes about two-thirds of the Cen peak excess, while the peak emission from the Perseus–Pisces supercluster is almost negligible.

In summary, by simply requiring the presence of hot spots in the UMa and Cen regions, the model accurately reproduces the patterns of flux over- and underdensities across the sky as seen by the Pierre Auger and Telescope Array observatories. The model in Figure 5 also naturally reproduces the amplitude of the flux excesses in the hot spot regions, as well as the angular correlation scale, $\Theta_{\text{obs}}$, for the magnetic-field parameters in the Local Sheet in Equation (13).

**4. Astrophysical Sources of Ultrahigh-energy Cosmic Rays**

**4.1. *Transient Rather Than Persistent Sources***

We have hypothesized that the sources responsible for the
acceleration of UHECRs and for the localized excesses in the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Model flux map of UHECRs at energies above 40 EeV for a burst rate per unit SFR $k = 1 \times 10^{-4} M_{\odot}^{-1}$. The sky map is smoothed using a top-hat function of radius $\Psi = 25^{\circ}$. The red line shows the supergalactic plane. The area delimited by the gray line represents the zone of avoidance. Following Biteau (2021), prominent galaxies and structures contributing to the flux are represented by gray markers and the center of the Local Void is indicated by a white circle. (The complete figure set (seven images) is available in the online article.)</figcaption></figure>

sky are transient. This hypothesis stems from a criterion by Hillas (1984), which requires UHECRs to remain confined for at least one gyration within the emitter. The Hillas criterion can be neatly reformulated by noting that the vacuum impedance $\Omega_{vac} = (\epsilon_0 c)^{-1} = 377 \Omega$ imposes losses at the source of the order of (Lovelace 1976; Waxman 1995; Blandford 2000)

$$L_B > \frac{\Gamma^2 R^2}{\beta \Omega_{vac}} > 6.6 \cdot 10^{43} \text{ erg s}^{-1} \times \left(\frac{\Gamma^2/\beta}{100}\right) \left(\frac{R}{5 \text{ EV}}\right)^2, \quad (17)$$

where $L_B$ is the bolometric luminosity and $(\Gamma, \Gamma\beta)$ is the four velocity of the ejecta within which acceleration takes place. Maintaining such a nonthermal luminosity over a long period of time is difficult. For example, the magnetized winds of the starburst galaxy M82 carry only $5 \times 10^{41} \text{ erg s}^{-1}$ (Contursi et al. 2013). Even the bolometric emission of M82, $2.5 \times 10^{44} \text{ erg s}^{-1}$ essentially of thermal origin, only slightly exceeds the Lovelace–Waxman–Blandford limit in Equation (17).

The evidence of anisotropy correlated with galaxies in the Local Sheet could come instead from the echo of a past burst from the jetted active galactic nucleus of Cen A (deboosted nonthermal bolometric luminosity of $10^{43} \text{ erg s}^{-1}$; Borkar et al. 2021). The excesses of UHECRs in the sky would then be due to reflections of the past burst off the magnetized winds of M82 and off the other starburst galaxies in the Local Sheet (Bell & Matthews 2022). For the viability of this scenario, heavy nuclei must be injected into the jets of Cen A and other active galactic nuclei (see, e.g., the “espresso” model of Mbarek & Caprioli 2021). Whether such scenarios meet the production rate requirements remains to be determined. The alternative that can be investigated in a more quantitative way is a transient

origin of UHECRs from stellar explosions within all star-forming galaxies, as discussed here.

## 4.2. Stellar Rather Than Galactic Size

To evaluate the viability of any astrophysical scenario involving ultrahigh-energy nuclear emission, the metallicity of the accelerated ejecta must be examined. The primary or secondary origin of the protons detected around the ankle of the UHECR energy spectrum depends on the evolution of the injection with redshift (see Section 3.1). Thus, the metallicity of the source ejecta, i.e., the ratio of nuclei heavier than helium to protons, is poorly constrained. The ratios of helium to metals and of heavy to intermediate-mass nuclei prove to be more constraining observables.

Because of the very hard spectrum of the nuclei, of index $\gamma_A < 1$, the integral of the emitted spectrum is convergent: the number of emitted nuclei of each species is proportional to $\epsilon_A/E_{\text{max}}^{Z_A}$. The mass fractions of He to metals is estimated from the best-fit parameters provided in Section 3.1 as

$$\begin{aligned} \frac{M(\text{He})}{M(\text{C-Fe})} \bigg|_{\text{UHECR}} &= \frac{A_{\text{He}} \frac{\epsilon_{\text{He}}}{E_{\text{max}}^{Z_{\text{He}}}}}{A_{\text{CNO}} \frac{\epsilon_{\text{CNO}}}{E_{\text{max}}^{Z_{\text{CNO}}}} + A_{\text{Si}} \frac{\epsilon_{\text{Si}}}{E_{\text{max}}^{Z_{\text{Si}}}} + A_{\text{Fe}} \frac{\epsilon_{\text{Fe}}}{E_{\text{max}}^{Z_{\text{Fe}}}}} \\ &= \frac{\epsilon_{\text{He}}}{\epsilon_{\text{CNO}} + \epsilon_{\text{Si}} + \epsilon_{\text{Fe}} \frac{A_{\text{Fe}}}{2Z_{\text{Fe}}}} \\ &= 0.21 \pm 0.05_{\text{stat.}} \pm 0.06_{\text{sys.}}, \end{aligned} \quad (18)$$

where the first uncertainty, being statistical, takes into account the covariance between the nuclear components and where the second, being systematic, characterizes the variations in energy scale, hadronic-interaction generator, and redshift evolution. The value inferred for UHECRs is 2 orders of magnitude below
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** UHECR burst rate as a function of the cosmic-ray energy emitted per burst (left), and as a function of the kinetic energy of the outflow (right). The solid red and dotted gold lines delimit the regions allowed by the composition, spectrum, and arrival-direction data when including or not transient emission from the Milky Way in the model. On the left-hand side, the colors of the two y-axes show the burst rate per unit matter obtained for UHECR sources evolving as SFRD or SMD. On the right-hand side, the allowed regions are based on an efficiency of conversion of kinetic energy to particles, $\eta$, between 1% and 10%. Rectangular markers indicate the statistical uncertainty associated with the number of X-ray bursts observed for each type of source. The vertical error bars show the range associated with the beaming correction factor, while the horizontal error bars illustrate the range of kinetic energy per burst in each population. The color bar shows the threshold luminosity above which each source density is measured. The low luminosities excluded by the Hillas–Lovelace–Waxman–Blandford criterion is shown as a green hatched region of the color bar.</figcaption></figure>

<table><caption>**Table 2**

Properties of the Stellar-sized X-Ray Transients</caption><thead><tr><th>Type</th><th>Beam Angle (deg)</th><th>log<sub>10</sub>E<sub>K</sub> (erg)</th><th>References</th><th>log<sub>10</sub>L<sub>min</sub> (erg s<sup>-1</sup>)</th><th>log<sub>10</sub>n&#x0307;<sub>obs</sub> (Gpc<sup>-3</sup> yr<sup>-1</sup>)</th><th>References</th><th>log<sub>10</sub>n&#x0307;<sub>true</sub> (Gpc<sup>-3</sup> yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>l/L-GRBs</td><td>5–20</td><td>48–50</td><td>Cano et al. (2017)</td><td>46</td><td>2.64 ± 0.21</td><td>Sun et al. (2015)</td><td>4.5 ± 0.6</td></tr><tr><td>SBOs</td><td>...</td><td>50–52</td><td>Waxman & Katz (2017)</td><td>44</td><td>4.24 ± 0.39</td><td>Sun et al. (2015)</td><td>4.2 ± 0.4</td></tr><tr><td>hL-GRBs</td><td>1–5</td><td>50–53</td><td>Cano et al. (2017)</td><td>50</td><td>-0.10 ± 0.06</td><td>Sun et al. (2015)</td><td>3.0 ± 0.7</td></tr><tr><td>TDEs</td><td>...</td><td>48–50</td><td>Cendes et al. (2022)</td><td>44</td><td>2.67 ± 0.24</td><td>Sun et al. (2015)</td><td>2.7 ± 0.2</td></tr><tr><td>S-GRBs</td><td>5–20</td><td>50–52</td><td>Laskar et al. (2022)</td><td>50</td><td>0.12 ± 0.12</td><td>Sun et al. (2015)</td><td>1.9 ± 0.6</td></tr><tr><td>UL-GRBs</td><td>5–20</td><td>52–53</td><td>Beniamini et al. (2015)</td><td>48</td><td>-0.61 ± 0.39</td><td>Prajs et al. (2017)</td><td>1.2 ± 0.7</td></tr><tr><td>Jetted TDEs</td><td>5–20</td><td>51–52</td><td>Cendes et al. (2022)</td><td>48</td><td>-1.58 ± 0.42</td><td>Sun et al. (2015)</td><td>0.2 ± 0.7</td></tr></tbody></table>

**Note.** In column (8), the true rate density of bursts $\dot{n}_{true} = f_b^{-1}\dot{n}_{obs}$ is determined from the rate observed in soft and hard X-rays, $\dot{n}_{obs}$ (column (6)), above a luminosity threshold $L_{min}$ (column (5)). The true rate density also accounts for the beaming correction factor of the relativistic component, $f_b = 1 - \cos\theta_{jet}$, where $\theta_{jet}$ is the two-sided jet opening angle, or beam angle (column (2)). The latter is taken as 90° for noncollimated outflows from SBOs and TDEs.

$M(\text{He})/M(\text{C-Fe})|_{\text{ISM}} = 18 \pm 2$ (Lodders et al. 2009), which rules out an ejecta drawn from the interstellar medium.

The mass fraction of heavy nuclei (above Ne) to intermediate nuclei (He and CNO) is estimated as

$$ \frac{M(\text{Ne-Fe})}{M(\text{He-O})} \bigg|_{\text{UHECR}} = \frac{\varepsilon_{\text{Si}} + \varepsilon_{\text{Fe}} \frac{A_{\text{Fe}}}{2Z_{\text{Fe}}}}{\varepsilon_{\text{He}} + \varepsilon_{\text{CNO}}} \\ = 0.30 \pm 0.05_{\text{stat.}} \pm 0.10_{\text{sys.}} \quad (19) $$

A stellar-sized transient scenario should then account for a proportion of Ne–Fe elements (synthesized by C, O, and Si burning) to CNO elements (synthesized by the CNO cycle and He burning) at least as large as the value in Equation (19). This is in line with the ratio of heavy to CNO nuclei dispersed by the explosions of massive stars in the interstellar medium, namely $M(\text{Ne-Fe})/M(\text{CNO}) = 0.53 \pm 0.09$ (Lodders et al. 2009; see also Zhang et al. 2018 for more detailed accretion-ejection models).<fnref n="8" />

The elemental distribution at the sources appears to be consistent with that of massive stars stripped of their outer H

and He layers. On the other hand, the He-poor composition inferred in Equation (18) rules out the acceleration of ionized atoms from the interstellar or intergalactic medium, as might be expected, e.g., in the termination shocks of jets from active galactic nuclei (see, e.g., Matthews et al. 2019). The composition argument alone does not exclude the acceleration, within these jets, of ejecta from the disruption or explosion of a massive star (see, e.g., Zhang et al. 2017; Biehl et al. 2018b; Bosch-Ramon 2023, for recent references). Consideration of the rate and energetics of UHECR bursts provided in the next section, however, narrows down the populations of astrophysical sources that are candidates for ultrahigh-energy acceleration, at least for those populations whose burst-rate density is well constrained by observations.

### 4.3. Burst Rate and Kinetic Energy

As shown in Figure 4, arrival-direction data constrain the UHECR burst rate per unit SFR or $M_\star$. Such constraints are illustrated by the horizontal bands in Figure 6, left, which include systematic uncertainties in the SFR and $M_\star$. The constraints are fairly independent of the tracer normalized relative to the SFR and $M_\star$ density measurements by Driver et al. (2018) up to a redshift $z=0.08$. A more direct comparison between the model and the UHECR arrival

<footnote n="8">⁸ The systematic uncertainty on the quantity estimated in Equation (19) is smaller than that on the $M(\text{Ne-Fe})/M(\text{CNO})$ ratio, as can be seen by comparing the CNO proportions for the SMD and SFRD evolution scenarios.</footnote>
directions than that presented in Section 3.3, e.g., using an approach similar to that employed by Aab et al. (2018) and Abreu et al. (2022), is expected not only to yield a correlated signal fraction larger than that inferred with the more limited catalogs tested thus far (all the galaxies in the volume are virtually considered in the approach presented here) but also to restrict even more tightly the parameter space. The latter is limited here by the mere existence of hot spots in the two prominent sky regions identified by the Pierre Auger and Telescope Array collaborations.

In Figure 6, left, the diagonal bands show the constraints on the energy production rate, $\mathcal{L} \propto k\mathcal{E}_{\text{UHECR}}$ with $\mathcal{E}_{\text{UHECR}} = \sum_A \mathcal{E}_A$, as derived from the fit to the spectral and composition data (see Section 3.1). These diagonal bands take into account the systematic uncertainty in the energy scale, which is a dominant factor due to the rapid decrease of the UHECR flux with energy. Combining the constraints from the spectral and composition data with those from arrival directions partially breaks the degeneracy between burst-rate density and the quantity of UHECRs emitted per burst. The resulting constraints are represented by the gold and red bands in Figure 6, left, which are obtained by including or not including transient emission from the Milky Way in the UHECR sky model, respectively (see Section 3.3).

Assuming a conversion efficiency of the kinetic energy or Poynting flux of the ejecta into cosmic rays, $\eta = 1\%-10\%$, we bracket in Figure 6, right, the burst-rate density and kinetic energy per burst to constrain the populations of candidate sources. We use the kinetic energy, $E_K$, of the ejecta rather than its electromagnetic fluence, since the former is more directly related to the total cosmic-ray energy of each burst through $\mathcal{E}_{\text{UHECR}} = \eta E_K$.

Soft and hard X-ray observations, mainly from Swift, have been instrumental in characterizing the properties of the populations of stellar-sized explosions and mergers listed in Table 2, namely supernova shock breakouts (SBOs), nonjetted and jetted tidal disruption events (TDEs), short-duration gamma-ray bursts (S-GRBs), and L-GRBs. The latter are divided in Figure 6 into subcategories: high-luminosity and low-luminosity (lL-GRBs and hL-GRBs) for L-GRBs above and below $10^{50}$ erg s$^{-1}$, respectively, as well as ultralong GRBs (UL-GRBs), whose prompt emission lasts several kiloseconds (see references in Table 2). As indicated by the color bar in Figure 6, right, the luminosity of nonjetted TDEs and SBOs hardly meets the Hillas–Lovelace–Waxman–Blandford criterion in Equation (17). Jetted TDEs appear to be too infrequent to explain the UHECR sky. The kinetic energy of the ejecta of S-GRBs appears to be too low, although tighter constraints on their outflows may be needed to formally settle the matter. The various classes of L-GRBs seem to satisfy all the criteria presented in Figure 6, right. We note that hL-GRBs occupy a central position in the parameter space.

## 5. Conclusion

In this study, we develop a three-dimensional model of the UHECR sky aimed at reproducing the main features of the measurements from the Pierre Auger Observatory and Telescope Array in the flux-suppression region, i.e., above 40 EeV. In this energy range, excesses on Gaussian angular scales of the order of $15^\circ$ are emerging, as indicated by pieces of evidence between $4.2\sigma$ and $4.6\sigma$ provided by the Pierre Auger Collaboration (Abreu et al. 2022) and the joint working group with the Telescope Array Collaboration (Abdul Halim et al. 2023).

The evidence comes from a 10%-20% anisotropic excess of events correlated with a flux-limited catalog of about 40 star-forming galaxies. The model proposed here takes advantage of the 400,000-galaxy catalog of Biteau (2021) and of the tensor propagation formalism discussed in Section 3.1, which allow us to efficiently capture close to 100% of the UHECR flux expected above 40 EeV, irrespective of the cosmological evolution of the sources. The fit of a generic model to composition and spectral data from $\simeq 1$ to $\simeq 100$ EeV constrains the UHECR emissivity per unit matter as well as the ratio of heavy elements to intermediate-mass nuclei at escape from the sources. The latter suggests the ejecta composition of short-lived massive stellar progenitors, which supports an UHECR origin in transient stellar-sized sources. The software that enables us to efficiently fit spectral and composition data and to produce flux maps is made publicly available (Condorelli et al. 2024).

In addition to the three spatial dimensions covered with unprecedented depth at ultrahigh energies, this study also accounts for the temporal dimension. We consistently consider the temporal and angular spreads caused by the magnetic fields that UHECRs cross. We identify the parameter space region where the median transient sky model shows an excess in the Cen region and an excess in the UMa region, as suggested by UHECR observations above 40 EeV. The flux contrast and angular spread inferred from the observations are found in this region of parameter space.

The nodal point of the model, which makes it falsifiable, is the existence of a turbulent magnetic field in the Local Sheet. Its coherence length of 10 kpc and amplitude of 0.5–20 nG match the properties inferred for galaxy filaments (Vernstrom et al. 2021; Carretti et al. 2022). Whether radio observatories such as LOFAR could already probe the signatures of such a field, or whether future measurements with the Square Kilometer Array could reveal it, remains to be determined (see, e.g., Akahori et al. 2018 for a review). Our elementary model deliberately leaves out the coherent fields, the specific realization of which necessarily affects the comparison between model and UHECR data (e.g., Bray & Scaife 2018; Hackstein et al. 2018; Ding et al. 2021; Unger & Farrar 2024, for recent discussions). We discuss an example of factoring in coherent Galactic deflection models in the Appendix.

The presence of excesses in the Cen and UMa regions allows us to bracket the burst rate per unit matter. This, combined with the production rate of UHECRs, constrains the total energy emitted by each burst of UHECRs and hence the kinetic energy of the outflows within which acceleration takes place. The conclusion is that, among the stellar-sized transient X-ray sources with a measured rate density, only the L-GRB population can be the source of all UHECRs. The discovery of the sources requires the confirmation of the anisotropy of the UHECR sky above 40 EeV with $5\sigma$ confidence level.

More precise measurements of the properties of populations of nonthermal transient sources will certainly yield a sharper picture of the L-GRB subtype favored for UHECR acceleration. Unprecedented exploration of transient populations is expected with the development of time-domain astronomy, both at the highest observed photon energies with the Cherenkov Telescope Array Observatory, in X-rays with satellites such as Einstein Probe, and at optical wavelengths with the Vera C. Rubin Observatory. These observations will also shed light on the acceleration, radiation, and escape
processes that take place in transient astrophysical sources (see, e.g., Murase & Bartos 2019, and references therein), subjects we have deliberately left aside in this macroscopic study of the cosmic-ray sky at the highest energies.

We have limited the scope of this work to modeling the cosmic-ray sky in the flux-suppression region beyond 40 EeV, where our galaxy catalog captures over 97% of the expected flux. The lower-energy region, down to the ankle of the spectrum at 5 EeV, also presents particularly interesting observational features (Aab et al. 2017b, 2020). The observation of a flux dipole of increasing amplitude from 4–8 EeV to 16–32 EeV has been the subject of numerous models, all suggesting qualitative agreement with a distribution of sources following that of large-scale structures in the local Universe (Harari et al. 2015; Hackstein et al. 2018; Wittkowski & Kampert 2018; Allard et al. 2022; Bister & Farrar 2024). Our model, like the others mentioned above, yields an increase in dipole amplitude with increasing energy. However, measurements suggest that cosmic-ray rigidity is 2–3 times lower at 8 EeV than in the flux-suppressed region (Alves Batista et al. 2019), so that angular deviations are expected to be 2–3 times larger. The coherent components of the magnetic fields traversed therefore appear to be an essential ingredient for fully modeling the UHECR dipole. In addition to the advances on this front from radio astronomy, the SPHEREx infrared satellite has the potential to provide a complete spectroscopic mapping of the nearby Universe, and thus to deepen the reconstruction of the source distribution responsible for the UHECR dipole.

The constraints drawn on extragalactic UHECR accelerators may also have consequences for cosmic-ray nuclei observed below the ankle energy. With rare transient events in each galaxy, we have seen that the nonobservation of trans-EeV particles from past acceleration episodes in the Milky Way further shrinks the allowed range of burst rate per unit matter. As the SFR of the Milky Way is estimated at $SFR_{MW} \simeq 1.6 M_{\odot} \text{yr}^{-1}$ (Licquia & Newman 2015), our constraints result in a Galactic burst rate of $kSFR_{MW} \simeq 1 \times 10^{-5}-3 \times 10^{-3} \text{yr}^{-1}$, which brackets the estimates reported in Loeb & Waxman (2002) and Katz et al. (2013). Further investigations of the Galactic confinement of cosmic rays are needed to explore the viability of rare Galactic events to fuel the energy density of subankle nuclei; they are left for future studies.

## Acknowledgments

We acknowledge the work of the last two reviewers. Their constructive comments helped to improve the form and content of this manuscript. We are also very grateful to the colleagues who kindly provided comments on this manuscript, in particular Martin Lemoine, Kohta Murase, Peter Tyniakov, and Roger Clay, as a member of the Pierre Auger Collaboration publication committee. We gratefully acknowledge funding from ANR via the grant Multi-messenger probe of Cosmic Ray Origins (MICRO), ANR-20-CE92-0052. This work was also made possible with the support of the Institut Pascal at Université Paris-Saclay during the Paris-Saclay Astroparticle Symposium 2022, with the support of the P2IO Laboratory of Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of Université Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3

master projet UCMN, EuCAPT ANR-11-IDEX-0003-01, and Paris-Saclay ANR-10-LABX-0038.

# Appendix
## Impact of the Galactic Magnetic Field

When constructing the sky maps, the impact of the magnetic fields has been modeled by spreading the flux from a source over an angular scale $\Theta$. This is an approximation that neglects mass-spectrometric-type as well as flux-magnification effects, which can become important particularly in the case of large-scale fields such as the Galactic one. We show here the extent to which this approximation proves relevant.

We follow Jansson & Farrar (2012) for the Galactic magnetic field and consider Kolmogorov turbulence with an intensity scaled to the coherent component, and turbulence alone for the Local Sheet. The deflections from the turbulence follow from Achterberg et al. (1998). To build the sky maps, we perform numerical experiments using $10^6$ test particles. Cosmic-ray particles with opposite charge are backtracked from their arrival direction $n_k$ observed on Earth, drawn at random from a uniform distribution on the sphere. Arrival directions $n_k$ are then given a weight proportional to those of the sources sampled within $0.5^\circ$ along the trajectory of the test particles. Using a backtracking approach, the directional-dependent effects of flux magnification, defined as the flux received from a given source relative to that received in the absence of the magnetic field, are handled automatically. In this respect, our results are consistent with those of Farrar (2014) and Farrar & Sutherland (2019). This would not be the case should we forward-track particles from sources to Earth by means of lookup tables to map the entry point in the Galaxy with the arrival direction of a single UHECR on Earth, unless magnification factors are introduced into the lookup tables as in Bretz et al. (2014). We note that the magnification effects stemming from the turbulence may be overlooked here. The geometry of the turbulence in the Galaxy may in addition not be accurately described by a homogeneous and isotropic power spectrum, as an anisotropic one (Goldreich & Sridhar 1995) has been shown to describe better, at least locally, the cosmic-ray anisotropies in the TeV–PeV energy range (Giacinti & Kirk 2017). A full account of these effects requires to track test particles in specific realizations of the turbulence, beyond the application of random deflections. Such a comprehensive study is beyond the scope of the present work and will be reported elsewhere.

The exact content in each element received on Earth impacts the results due to the rigidity dependencies of both the coherent deflections and of the flux-magnification effects in a large-scale structured magnetic field. We consider for simplicity the benchmark mass composition as inferred from the best-fit parameters. We show the results of this procedure in Figure 7. A comparison with Figure 5 shows that the main changes induced by the large-scale component of the Galactic magnetic field are a flux magnification of the UMa region (around M81/82) and a pattern distortion along the supergalactic plane from NGC 1569 to NGC 253. With caution regarding both the model of the Galactic magnetic field and the flux determination currently obtained in the UMa region with the Telescope Array exposure, we note that flux magnification from the model results in a flux contrast between the UMa and Cen regions that may better match the observations. On the other hand, for a burst rate $k = 1 \times 10^{-4} M_{\odot}^{-1}$, the model suggests a flux
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Model flux map of UHECRs at energies above 40 EeV including deflections in the model of Jansson & Farrar (2012) of the Galactic magnetic field, for a burst rate per unit SFR $k = 1 \times 10^{-4} M_{\odot}^{-1}$. The sky map is smoothed with a top-hat function of radius $\Psi = 25^{\circ}$ to ensure a meaningful comparison with the observations. The red line shows the supergalactic plane. The area delimited by the gray line represents the zone of avoidance.</figcaption></figure>

excess between NGC 1569 and NGC 253 at the level of $\simeq 12 \times 10^{-3} \text{ km}^{-2} \text{ yr}^{-1} \text{ sr}^{-1}$, which remains difficult to uncover with the current exposures of the Pierre Auger Observatory and Telescope Array in this part of the sky.

The exploratory study case discussed here shows the expected impact of the large-scale component of the Galactic magnetic field, the knowledge of which still suffers from large uncertainties. It is observed that both mass-spectrometric-type and flux-magnification effects yield a sky map that is consistent with the current observational data.

## ORCID iDs

Sullivan Marafico https://orcid.org/0000-0003-4604-4677

Jonathan Biteau https://orcid.org/0000-0002-4202-8939

Antonio Condorelli https://orcid.org/0000-0001-5681-0086

Olivier Deligny https://orcid.org/0000-0001-6863-6572

Johan Bregeon https://orcid.org/0000-0002-6790-5328

## References

* Aab, A., Abreu, P., Aglietta, M., et al. 2014a, PhRvD, 90, 122005

* Aab, A., Abreu, P., Aglietta, M., et al. 2014b, PhRvD, 90, 122006

* Aab, A., Abreu, P., Aglietta, M., et al. 2015, NIMPA, 798, 172

* Aab, A., Abreu, P., Aglietta, M., et al. 2017a, JCAP, 04, 038

* Aab, A., Abreu, P., Aglietta, M., et al. 2017b, Sci, 357, 1266

* Aab, A., Abreu, P., Aglietta, M., et al. 2018, ApJL, 853, L29

* Aab, A., Abreu, P., Aglietta, M., et al. 2020, ApJ, 891, 142

* Abdul Halim, A., Abre, P., Aglietta, M., et al. 2023, ICRC (Nagoya), 444, 521

* Abdul Halim, A., Abreu, P., Aglietta, M., et al. 2024, JCAP, 2024, 022

* Abreu, P., Aglietta, M., Albury, J. M., et al. 2021, EPJC, 81, 966

* Abreu, P., Aglietta, M., Albury, J. M., et al. 2022, ApJ, 935, 170

* Abu-Zayyad, T., Aida, R., Allen, M., et al. 2013, NIMPA, 689, 87

* Achterberg, A., Gallant, Y. A., Norman, C. A., & Melrose, D. B. 1998, in 19th Texas Symp. Relativistic Astrophysics and Cosmology, ed. J. Paul, T. Montmerle, & E. Aubourg, 617

* Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32

* Aharonian, F., Aschersleben, J., Backes, M., et al. 2023, ApJL, 950, L16

* Akahori, T., Nakanishi, H., Sofue, Y., et al. 2018, PASJ, 70, R2

* Allard, D., Aublin, J., Baret, B., & Parizot, E. 2022, A&A, 664, A120

* Allard, D., Aublin, J., Baret, B., & Parizot, E. 2024, A&A, 686, A292

* Aloisio, R., Berezinsky, V., & Blasi, P. 2014, JCAP, 10, 020

* Aloisio, R., Boncioli, D., Di Matteo, A., et al. 2017, JCAP, 11, 009

* Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, FrASS, 6, 23

* Aragon-Calvo, M. A., Silk, J., & Neyrinck, M. 2023, MNRAS, 520, L28

* Bell, A. R., & Matthews, J. H. 2022, MNRAS, 511, 448

* Beniamini, P., Nava, L., Duran, R. B., & Piran, T. 2015, MNRAS, 454, 1073

* Biehl, D., Boncioli, D., Fedynitch, A., & Winter, W. 2018a, A&A, 611, A101

* Biehl, D., Boncioli, D., Lunardini, C., & Winter, W. 2018b, NatSR, 8, 10828

* Bister, T., & Farrar, G. R. 2024, ApJ, 966, 71

* Biteau, J. 2021, ApJS, 256, 15

* Blandford, R. D. 2000, PhST, 85, 191

* Bonafede, A., Feretti, L., Murgia, M., et al. 2010, A&A, 513, A30

* Borkar, A., Adhikari, T. P., Różańska, A., et al. 2021, MNRAS, 500, 3536

* Bosch-Ramon, V. 2023, A&A, 677, L14

* Bray, J. D., & Scaife, A. M. M. 2018, ApJ, 861, 3

* Bretz, H.-P., Erdmann, M., Schiffer, P., Walz, D., & Winchen, T. 2014, APh, 54, 110

* Cano, Z., Wang, S.-Q., Dai, Z.-G., & Wu, X.-F. 2017, AdAst, 2017, 8929054

* Carretti, E., Vacca, V., O'Sullivan, S. P., et al. 2022, MNRAS, 512, 945

* Cendes, Y., Berger, E., Alexander, K. D., et al. 2022, ApJ, 938, 28

* Condorelli, A., Biteau, J., & Adam, R. 2023, ApJ, 957, 80

* Condorelli, A., Marafico, S., Biteau, J., & Bregeon, J. 2024, Software and dataset for the paper "Closing the Net on Transient Sources of Ultra-high Energy Cosmic Rays," v3.0, Zenodo, doi:10.5281/zenodo.11440864

* Contursi, A., Poglitsch, A., Graciá Carpio, J., et al. 2013, A&A, 549, A118

* di Matteo, A., Bister, T., Biteau, J., et al. 2020, ICRC (Madison, WI), 36, 439

* Ding, C., Globus, N., & Farrar, G. R. 2021, ApJL, 913, L13

* Dolag, K., Kachelrieb, M., & Semikoz, D. V. 2009, JCAP, 01, 033

* Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891

* Ehlert, D., Oikonomou, F., & Unger, M. 2023, PhRvD, 107, 103045

* Erdmann, M., Müller, G., Urban, M., & Wirtz, M. 2016, APh, 85, 54

* Fang, K., & Murase, K. 2018, NatPh, 14, 396

* Farrar, G. R. 2014, CRPhy, 15, 339

* Farrar, G. R., & Sutherland, M. S. 2019, JCAP, 05, 004

* Fitoussi, T., Abdul Halim, A., Abreu, P., et al. 2023, ICRC (Nagoya), 444, 319

* Giacinti, G., & Kirk, J. G. 2017, ApJ, 835, 258

* Gilmore, R. C., Somerville, R. S., Primack, J. R., & Dominguez, A. 2012, MNRAS, 422, 3189

* Globus, N., Allard, D., & Parizot, E. 2015, PhRvD, 92, 021302

* Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

* Greisen, K. 1966, PhRvL, 16, 748

* Hackstein, S., Vazza, F., Brüggen, M., Sorce, J. G., & Gottlöber, S. 2018, MNRAS, 475, 2519
Harari, D., Mollerach, S., & Roulet, E. 2015, PhRvD, 92, 063014

Harari, D., Mollerach, S., & Roulet, E. 2016, JCAP, 2016, 010

Hillas, A. M. 1984, ARA&A, 22, 425

Hoffman, Y., Pomarede, D., Brent Tully, R., & Courtois, H. 2017, NatAs, 0036

Huchra, J. P., Macri, L. M., Masters, K. L., et al. 2012, ApJS, 199, 26

Jansson, R., & Farrar, G. R. 2012, ApJ, 757, 14

Jedamzik, K., & Saveliev, A. 2019, PhRvL, 123, 021301

Kaiser, N. 1986, MNRAS, 222, 323

Katz, B., Waxman, E., Thompson, T., & Loeb, A. 2013, arXiv:1311.0287

Koning, A. J., Hilaire, S., & Duijvestijn, M. C. 2005, in AIP Conf. Proc. 769, 1154

Koning, A. J., & Rochman, D. 2012, NDS, 113, 2841

Kotera, K., Allard, D., Murase, K., et al. 2009, ApJ, 707, 370

Koushan, S., Driver, S. P., Bellstedt, S., et al. 2021, MNRAS, 503, 2033

Laskar, T., Escorial, A. R., Schroeder, G., et al. 2022, ApJL, 935, L11

Lavaux, G., & Hudson, M. J. 2011, MNRAS, 416, 2840

Lemoine, M., Sigl, G., Olinto, A. V., & Schramm, D. N. 1997, ApJL, 486, L115

Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96

Lodders, K., Palme, H., & Gail, H.-P. 2009, LanB, 4B, 712

Loeb, A., & Waxman, E. 2002, arXiv:astro-ph/0205272

López Fernández, R., González Delgado, R. M., Pérez, E., et al. 2018, A&A, 615, A27

Lovelace, R. V. E. 1976, Natur, 262, 649

Luce, Q., Marafico, S., Biteau, J., Condorelli, A., & Deligny, O. 2022, ApJ, 936, 62

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2019, MNRAS, 482, 4303

Mbarek, R., & Caprioli, D. 2021, ApJ, 921, 85

McCall, M. L. 2014, MNRAS, 440, 405

Miralda-Escude, J., & Waxman, E. 1996, ApJL, 462, L59

Murase, K., & Bartos, I. 2019, ARNPS, 69, 477

Murase, K., & Takami, H. 2009, ApJL, 690, L14

Neronov, A., & Vovk, I. 2010, Sci, 328, 73

O'Sullivan, S. P., Brüggen, M., Vazza, F., et al. 2020, MNRAS, 495, 2607

Peters, B. 1961, NCim, 22, 800

Pierog, T., Karpenko, I., Katzy, J. M., Yatsenko, E., & Werner, K. 2015, PhRvC, 92

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 596, A101

Prajs, S., Sullivan, M., Smith, M., et al. 2017, MNRAS, 464, 3568

Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, PhRvL, 116, 191302

Riehn, F., Engel, R., Fedynitch, A., Gaisser, T. K., & Stanev, T. 2020, PhRvD, 102, 063002

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40

Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJL, 747, L31

Sigl, G. 2001, in Physics and Astrophysics of Ultra-High-Energy Cosmic Rays, ed. M. Lemoine & G. Sigl (Berlin Heidelberg: Springer), 196

Sigl, G., Lemoine, M., & Olinto, A. V. 1997, PhRvD, 56, 4470

Sironi, L., Keshet, U., & Lemoine, M. 2015, SSRv, 191, 519

Sun, H., Zhang, B., & Li, Z. 2015, ApJ, 812, 33

Supanitsky, A. D., Cobos, A., & Etchegoyen, A. 2018, PhRvD, 98, 103016

Tavecchio, F., Ghisellini, G., Foschini, L., et al. 2010, MNRAS, 406, L70

Taylor, A. M., Ahlers, M., & Hooper, D. 2015, PhRvD, 92, 063011

Tkachenko, O., Abdul Halim, A., Abreu, P., et al. 2023, ICRC (Nagoya), 444, 438

Unger, M., & Farrar, G. 2024, ApJ, 970, 95

Unger, M., Farrar, G. R., & Anchordoqui, L. A. 2015, PhRvD, 92, 123001

Urban, O., Simionescu, A., Werner, N., et al. 2014, MNRAS, 437, 3939

Vazza, F., Brüggen, M., Gheller, C., et al. 2017, CQG, 34, 234001

Vernstrom, T., Heald, G., Vazza, F., et al. 2021, MNRAS, 505, 4178

Watson, A. A. 2022, JHEAp, 33, 14

Waxman, E. 1995, PhRvL, 75, 386

Waxman, E., & Katz, B. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 967

Wittkowski, D., & Kampert, K.-H. 2018, ApJL, 854, L3

Zatsepin, G. T., & Kuz'min, V. A. 1966, JETPL, 4, 78

Zhang, B. T., Murase, K., Kimura, S. S., Horiuchi, S., & Mészáros, P. 2018, PhRvD, 97, 083010

Zhang, B. T., Murase, K., Oikonomou, F., & Li, Z. 2017, PhRvD, 96, 063007