# Ultra-high-energy Cosmic Ray Sources can be Gamma-ray Dim

Angelina Partenheimer, Ke Fang, Rafael Alves Batista, Rogerio Menezes de Almeida

► To cite this version:

Angelina Partenheimer, Ke Fang, Rafael Alves Batista, Rogerio Menezes de Almeida. Ultra-high-energy Cosmic Ray Sources can be Gamma-ray Dim. The Astrophysical Journal Letters, 2024, 967 (1), pp.L15. 10.3847/2041-8213/ad4359. hal-04563842

HAL Id: hal-04563842

https://hal.science/hal-04563842v1

Submitted on 27 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Ultra-high-energy Cosmic-Ray Sources Can Be Gamma-Ray Dim

Angelina Partenheimer¹ ID, Ke Fang¹ ID, Rafael Alves Batista²,³ ID, and Rogerio Menezes de Almeida⁴ ID

¹Department of Physics, Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI 53706, USA

²Instituto de Física Teórica UAM/CSIC, Calle Nicolás Cabrera 13-15, Cantoblanco, 28049 Madrid, Spain

³Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

⁴Instituto de Física, Universidade Federal do Rio de Janeiro, Caixa Postal 68528, Rio de Janeiro, RJ 21941-972, Brazil

Received 2024 February 23; revised 2024 April 2; accepted 2024 April 25; published 2024 May 17

## Abstract

Ultra-high-energy cosmic rays (UHECRs), accelerated hadrons that can exceed energies of $10^{20}$ eV, are the highest-energy particles ever observed. While the sources producing UHECRs are still unknown, the Pierre Auger Observatory has detected a large-scale dipole anisotropy in the arrival directions of cosmic rays above 8 EeV. In this work, we explore whether resolved gamma-ray sources can reproduce the Auger dipole. We use various Fermi Large Area Telescope catalogs as sources of cosmic rays in CRPropa simulations. We find that in all cases, the simulated dipole has an amplitude significantly larger than that measured by Auger, even when considering large extragalactic magnetic field strengths and optimistic source weighting schemes. Our result implies that the resolved gamma-ray sources are insufficient to account for the population of sources producing the highest-energy cosmic rays, and there must exist a population of UHECR sources that lack gamma-ray emission or are unresolved by the current-generation gamma-ray telescopes.

*Unified Astronomy Thesaurus concepts: Ultra-high-energy cosmic radiation (1733); Gamma-ray sources (633)*

## 1. Introduction

Ultra-high-energy cosmic rays (UHECRs), charged nuclei whose energies can exceed $10^{20}$ eV, are the highest-energy particles ever observed. The sources of these particles are still unidentified, and the mechanism by which these extremely energetic particles are accelerated remains a mystery (Kotera & Olinto 2011; Alves Batista et al. 2019). The Pierre Auger Observatory (Auger) has detected a dipole anisotropy in cosmic-ray arrival directions above 8 EeV. The measured dipole has an amplitude of 6.5% and points in a direction ~125° away from the Galactic center, supporting a hypothesis of extragalactic origin for the highest-energy cosmic rays. (Aab et al. 2017, 2018a, 2020a; Abdul Halim et al. 2023).

Previous cross-correlation studies have indicated some possible excesses of UHECRs around gamma-ray-producing sources. For example, the Pierre Auger Observatory recently detected a hotspot in the region of Centaurus A, NGC 4945, and M83 (Abreu et al. 2022). In addition, the Telescope Array hotspot for UHECRs > 57 EeV (Abbasi et al. 2014) in the northern sky lies in the region of a few nearby gamma-ray emitting sources, in particular Mkn 421, Mkn 180, and M82. Cross-correlation studies between gamma-ray-bright sources and UHECR events also found an excess of cosmic-ray events around starburst galaxies (Aab et al. 2018b).

Motivated by these observations, our work attempts to address whether the population of gamma-ray sources could encapsulate the population of UHECR-producing sources. We use extragalactic sources from the Fermi Large Area Telescope (LAT) gamma-ray source catalogs as the origins of cosmic rays in simulations with CRPropa3 (Alves Batista et al. 2016; Ajello et al. 2017; Abdollahi et al. 2020; Ajello et al. 2020; Alves Batista et al. 2022). We conclude that the known gamma-ray

sources are insufficient to reproduce a dipole as low as that observed by Auger, suggesting that there may exist a population of gamma-ray dim UHECR sources.

In Section 2 we discuss properties of gamma-ray catalogs with Fermi-LAT 4FGL-DR4 as a benchmark case, as well as preprocessing measures applied to the catalog before simulation. In Section 3 we describe simulations of UHECRs with CRPropa3, especially our approach to accounting for deflections in the extragalactic and Galactic magnetic fields. In Section 4 we present our simulation results, and in Section 5 we discuss the implication of our results.

## 2. Properties of the Fermi-LAT 4FGL-DR4 Catalog

The Large Area Telescope (LAT) on board the Fermi Satellite surveys the sky in the GeV energy range, and has produced several generations of gamma-ray source catalogs, including the LAT 14 yr Source Catalog (4FGL-DR4; Abdollahi et al. 2020; Ballet et al. 2023), Third Fermi-LAT Catalog of High-Energy Sources (3FHL; Ajello et al. 2017), and the Fourth Catalog of Active Galactic Nuclei (4LAC-DR3; Ajello et al. 2020). Below we use 4FGL-DR4, which is the latest catalog from Fermi-LAT in the 50 MeV–1 TeV energy range as the basis for our simulations. We present results using the other GeV gamma-ray catalogs, which are consistent with those using 4FGL-DR4, in Appendix C.

The 4FGL-DR4 catalog contains 7195 sources, some of which are Galactic. Motivated by the fact that the UHECR sources are likely extragalactic (Aab et al. 2017), we remove Galactic sources that are flagged as pulsars, supernova remnants, pulsar wind nebulae, globular clusters, binaries, star-forming regions, or novae. Many sources in the Fermi-LAT Galactic plane region are unassociated (Abdollahi et al. 2020), so in order to fully eliminate Galactic sources we must remove all unassociated low-latitude ($|b| < 10°$) sources. However, to verify that this cut does not significantly affect our simulated dipole, we perform an additional test where we leave a randomly selected set of low-latitude, unassociated sources chosen such that (1) the density per

<figure><img src="image_3.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.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sky map of the 4FGL-DR4 catalog after Galactic sources and unassociated sources within $|b| < 10^\circ$ are removed. The color bar indicates the gamma-ray energy flux in erg cm$^{-2}$ s$^{-1}$. The dimmest source in the catalog has energy flux $3.09 \times 10^{-13}$ erg cm$^{-2}$ s$^{-1}$.</figcaption></figure>

unit solid angle of low-latitude unassociated sources matches the density per solid angle of high-latitude unassociated sources and (2) the gamma-ray energy flux of all selected low-latitude unassociated sources falls within one standard deviation of the gamma-ray energy flux of the high-latitude unassociated sources. We confirm that there is no difference between this scenario and the one where all low-latitude unassociated sources are removed. In total, we remove 1740 sources from the catalog, 1146 of which are unassociated. A sky map of the resulting source catalog is shown in Figure 1.

The 4FGL catalog has no redshift data provided. However, the catalog shares sources with the 4LAC-DR3 (Ajello et al. 2020) and the 3FHL (Ajello et al. 2017) catalogs, which have redshift data provided for many of their sources. We obtain the redshift information for 1805 sources from 4LAC-DR3 and 24 additional (not in 4LAC-DR3) sources from 3FHL. It is noteworthy that this provides redshift data for some of the nearest and brightest 4LAC and 3FHL sources, which are significant contributors to the large-scale anisotropy that we evaluate for simulated UHECRs.

For those sources that do not have redshift data provided in 4LAC-DR3 or 3FHL, we assign random redshifts following the star formation rate (SFR) as described in Yüksel et al. (2008). Random redshifts are assigned up to a maximum redshift of $z=4$; this maximum value is chosen because it is the maximum redshift of sources in the Fermi-LAT 3FHL and 4LAC-DR3 catalogs. The source distribution before and after random redshift assignment is shown in Figure 2.

We note that the sources with redshift data provided are typically those that are the nearest and brightest, while sources without measured redshift tend to comprise a more dim, distant background. We quantify the uncertainty in our results due to fluctuations in random redshift assignment by performing simulations for several different realizations of the random background, and we verify that random redshift assignment does not significantly change our results. Finally, as a comparison to evaluate the impact of our choice to assign redshift following an SFR distribution, we test a scenario where sources without redshift data provided are assigned redshift following a uniform distribution, and we conclude that this does not meaningfully affect our simulated dipole.

### 3. Simulation

Without loss of generality, we consider two source weighting schemes. In the first scenario, we assume that the

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Redshift distribution of sources in 4FGL-DR4 after removing Galactic sources and low-latitude unassociated sources. Sources with known redshifts are indicated by the blue shaded region. Sources without distance information are assigned random redshifts following an SFR distribution (Yüksel et al. 2008) as indicated by the black dotted curve. The sum of the two are shown as the red solid curve.</figcaption></figure>

number of particles injected per source is weighted according to the gamma-ray energy flux at 100 MeV–100 GeV. We refer to this case as “gamma-ray flux weighting.” In the second case, we consider a scenario where the probability of detecting UHECRs at Earth coming from each source is the same. We call it “uniform weighting.” This scenario is likely unphysical as it assumes that more-distant sources are also more luminous. However, it can be justified because it accounts for the fact that at greater distances, the gamma-ray source catalogs are likely to be less complete. While neither of these weighting schemes are expected to be a perfect reflection of physical reality, considering these two very different scenarios is intended to provide perspective on the UHECR dipole attainable with the catalogs of resolved gamma-ray sources. While we ran preliminary tests using some other commonly considered weighting schemes, these weighting schemes yielded dipole amplitudes even larger than that found with the gamma-ray flux weighting, and were therefore redundant in the context of our results.

We simulate the propagation of cosmic rays from gamma-ray sources using CRPropa 3.2 (Alves Batista et al. 2016, 2022). We use an injection spectrum following a power law and exponential cutoff at maximum rigidity:

$$ \frac{dN}{dE} = \begin{cases} f_i E^{-\gamma}, & E \le ZR_{\max} \\ f_i E^{-\gamma} e^{1-E/ZR_{\max}}, & E > ZR_{\max} \end{cases} $$

We adopt a spectral index $\gamma = -1.82$ and maximum rigidity $R_{\max} = 10^{18.15}$ V. We adopt luminosity fractions $I_p = 0.620$, $I_{\text{He}} = 0.233$, $I_{\text{N}} = 0.533$, $I_{\text{Si}} = 0.148$, and $I_{\text{Fe}} = 0.024$ that fit the Auger unfolded spectrum reported in Aab et al. (2020b) and the measurements of the depth of shower maximum reported in Yushkov (2019) considering a scenario of sources distributed homogeneously and with the same luminosity. Using such a scenario results in a conservative estimate of the expected anisotropy, as a lighter composition or a higher rigidity cutoff would be expected to result in a stronger dipole (di Matteo & Tinyakov 2018).

Our simulation accounts for the energy losses due to interactions with the cosmic microwave background and extragalactic background light (Gilmore et al. 2012) via
electron pair production, photodisintegration, and photopion production, as well as nuclear decay.

We treat the extragalactic magnetic field (EGMF) as a random turbulent field with a Kolmogorov spectrum described by two parameters: rms value of the magnetic fields $B_{RMS}$, and maximum correlation length $l_{coh}$. In particular, we adopt $B_{RMS} \sqrt{l_{coh}} = 1 \times 10^{-13} \text{ G} \sqrt{\text{Mpc}}$. This choice is motivated by the upper bound on cosmological magnetic fields by Jedamzik & Saveliev (2019) ($B_{RMS} \lesssim 10^{-11} \text{ G}$). It is also in agreement with lower bounds derived from gamma-ray observations, which generally indicate $B_{RMS} \gtrsim 10^{-17} \text{ G}$ at coherence lengths $l_{coh} \gtrsim 100 \text{ kpc}$ (e.g., Fermi-LAT Collaboration 2018; MAGIC Collaboration 2023). Note that EGMF effects on UHECRs depend not only on $B_{RMS}$ but also on $l_{coh}$. However, limits on the coherence length are almost nonexistent or very weak (e.g., Alves Batista & Saveliev 2020).

The average deflection angle of a cosmic ray with rigidity R propagating a distance D in a magnetic field of rms strength B and coherence length $\lambda$ is given by Waxman & Miralda-Escude (1996):

$$ \theta_s \approx 0.025 \sqrt{\frac{D}{\lambda}} \left(\frac{\lambda}{10 \text{ Mpc}}\right) \left(\frac{B}{10^{-11} \text{ G}}\right) \left(\frac{R}{10^{20} \text{ V}}\right)^{-1}. \quad (1) $$

The angular distribution of cosmic rays after averaging over time is given by a simple Gaussian:

$$ P(\phi) = \frac{1}{\theta_s \sqrt{\pi}} e^{-\phi^2 / \theta_s^2}. \quad (2) $$

We simulate cosmic-ray propagation directly from source to observer, accounting for all interactions and energy losses that occur during propagation. We then assign the cosmic ray a deflection angle $\phi$ from the source following Equation (2) using its rigidity at injection. In addition to the deflection angle, the EGMF also introduces a time delay to particles traversing the field, $\tau \approx D\theta_s^2/(4c)$, corresponding to a distance $D_{add} = D\theta_s^2/4$. To account for this additional travel distance, each particle is set back a distance $D_{add}$ at its initial injection, allowing for any extra interactions and energy losses due to the time delay to occur.

To account for cosmic-ray deflections in the Galactic magnetic field (GMF), we adopt the GMF model of Jansson & Farrar (2012) and use the “lensing” technique from Alves Batista et al. (2016). We include both the random striated and turbulent components of the field model. For the turbulent field, we use Kolmogorov spectral index 5/3, coherence length 60 pc, minimum and maximum scales of the turbulence 8 and 272 pc, respectively, and field strength following the parameters in Table 1 of Jansson & Farrar (2012). To generate the lens, we perform backtracking simulations of antiprotons with discrete rigidities on $\log_{10}R/V \in [17.0, 21.0]$ and record their deflections in the GMF. These results are then combined to generate a Galactic lens. We apply the Galactic lens to simulated sky maps of UHECRs using CRPropa’s Parti-cleMapsContainer feature, which yields a normalized sky map of UHECR events.

The dipole amplitude and direction are evaluated using Rayleigh analysis following the method of the Pierre Auger Collaboration (Aab et al. 2017) and summarized in Appendix B.

<figure><img src="image_6.png" /><figcaption>Figure 3. Spectrum at Earth of simulated UHECRs from Fermi-LAT sources. The energy range below 8 EeV is shaded in gray as only events with energy greater than 8 EeV are used in evaluating the dipole as in Aab et al. (2017).</figcaption></figure>

<table><caption>Table 1
Dipole Amplitude and Direction for Simulated UHECRs from 4FGL-DR4 Sources in Both Weighting Scenarios</caption><thead><tr><th>Events &gt;8 EeV</th><th>Amplitude (%)</th><th>α<sub>d</sub> (deg)</th><th>δ<sub>d</sub> (deg)</th></tr></thead><tbody><tr><td>Uniform</td><td>33 ± 7</td><td>303 ± 17</td><td>-11 ± 9</td></tr><tr><td>Energy flux</td><td>74 ± 5</td><td>252 ± 7</td><td>62 ± 3</td></tr></tbody></table>

Note. The simulated data sets contain 32,000 events >8 EeV following Aab et al. (2017).

## 4. Results

We present our simulation results using the 4FGL-DR4 catalog in Section 4.1. We then consider the effect of unresolved gamma-ray sources in Section 4.2. Even though the uniform weight scenario is somewhat unphysical, we show the results from this case first. This is because, given the source distribution of the resolved gamma-ray sources, the uniform weighting scheme yields the best possible agreement with the dipole amplitude and composition measured by Auger. We show the results of the gamma-ray flux scenario, which has poorer agreement with the Auger observations, in Appendix A.

### 4.1. Simulated UHECRs from Fermi-LAT Sources

The simulated UHECR spectrum at Earth in the case of uniform source weighting is shown in Figure 3. The spectrum is in good agreement with the observed spectrum of Auger (Aab et al. 2020b). The spectrum in the case of gamma-ray flux weighting is shown in Appendix A.

The dipole amplitudes and corresponding directions from the simulated UHECRs are summarized in Table 1. A flux map corresponding to the uniform weighting case is shown in Figure 4; a flux map in the case of gamma-ray flux weighting is shown in Appendix A. The corresponding uncertainties are the standard deviations in the distribution of results from 10 separate simulations with different iterations of random redshift assignment.

In the case of gamma-ray flux weighting, we find a dipole amplitude of 74%, which is more than 10 times larger than that measured by Auger. In this case, the dipole is dominated by
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Flux map corresponding to the case of uniform source weighting. The flux map is obtained by taking the ratio of the smoothed count density and exposure map following Aab et al. (2017). The sky map is smoothed with a 45° top-hat function using the `healpy.sphtfnc.beam2bl` function.</figcaption></figure>

Markarian 421, an especially bright and relatively nearby source in the 4FGL-DR4 catalog. This disproportionately large effect by Markarian 421 is also noted in Abdul Halim et al. (2024). In the case of uniform source weighting, we evaluate a dipole amplitude of 33%. This weighting scheme allows the more distant sources in the catalog to generate a stronger background against the nearby sources. However, even in this case, the dipole amplitude is significantly larger than the 6.5% measured by Auger. In both scenarios, the dipole direction points to a different region of the sky from dipole direction measured by Auger.

In both scenarios, the very large dipole amplitude indicates that the density of resolved gamma-ray sources is too low to generate a dipole as low as that measured by Auger. This disagreement is further confirmed by the tension between the simulated and observed dipole directions. Our results imply that, given the magnetic field models used in this analysis, the resolved gamma-ray sources are insufficient to account for the population of UHECR sources producing the Auger dipole. In other words, there must be a population of UHECR sources that are unobserved by the current gamma-ray facilities, and are gamma-ray dim.

## 4.2. Unresolved Sources

As a final test, we consider the effect that unresolved gamma-ray sources could have on the simulated cosmic-ray dipole. We test unresolved source populations using two approaches. In the first approach, we infer an unresolved source population using a source count distribution and consider the UHECR dipole when these sources are weighted according to their energy flux. In the second approach, we generate a fixed number of unresolved sources that are weighted uniformly. In both cases, the unresolved sources are assigned positions following an isotropic distribution and random redshifts following the SFR distribution discussed previously. While it may be possible to distribute the unresolved source positions anisotropically such that the direction of the Auger dipole is reproduced, our assumption of an isotropic distribution is a simple, conservative approximation to yield the lowest possible dipole with the fewest number of additional sources. The unresolved sources are appended to the 4FGL-DR4 catalog,

and we evaluate the dipole of cosmic rays simulated from this new distribution.

For the case in which unresolved sources are weighted by energy flux, we base our unresolved source population on the source count distribution from Amerio et al. (2023). This work uses machine learning to derive a source count distribution following $\frac{dN}{dE} \propto S^{-2}$ below the 4FGL-DR4 sensitivity of $2 \times 10^{-10} \text{ ph cm}^{-2} \text{ s}^{-1}$. Following this distribution, we generate a population of unresolved sources extending to the minimum luminosity of $S_{\text{min}} = 5 \times 10^{-12} \text{ ph cm}^{-2} \text{ s}^{-1}$ used in the analysis of Amerio et al. (2023). We derive the energy flux of each source from its photon flux assuming a source energy spectrum $\frac{dN_e}{dE} \propto E^{-b_\gamma}$ with $b_\gamma = 2.23$, which is the median spectral index of 4FGL-DR4 extragalactic sources. We find that in this case, there is no measurable difference compared to the energy flux weighted simulation without the additional population of unresolved sources. This is due to the fact that, despite the very large number of added sources, the very low energy fluxes of the unresolved sources are not sufficient to balance out the effect of a bright source like Markarian 421.

In the second test, we incrementally add fixed numbers of unresolved sources to the 4FGL-DR4 catalog, and we weight all sources uniformly. In this case, we find that a number of ~80,000 added sources (corresponding to a local source density of ~$5 \times 10^{-4} \text{ Mpc}^{-3}$) is required to lower the dipole amplitude to a level analogous to that measured by Auger. This result is informative because it demonstrates that to have a dipole amplitude as low as that measured by Auger, the UHECR source density must be significantly larger than the resolved gamma-ray source population, especially the population of sources comprising the UHECR background. The result is also in agreement with other work (e.g., Bister & Farrar 2024) that has derived a similar source density necessary to achieve a dipole as low as Auger's.

The result of this test with added unresolved sources strengthens our conclusion that UHECR injection cannot be proportional to gamma-ray energy flux, and in fact that gamma-ray-dim sources must account for a substantial portion of UHECR accelerators. Assuming that the accepted EGMF and GMF models do not severely underestimate average magnetic field strengths, a dipole amplitude as low as Auger's can be achieved only with a very high source density, which the resolved gamma-ray sources are not sufficient to account for.

## 5. Summary and Discussion

In this work, we simulate UHECRs originating from extragalactic Fermi-LAT gamma-ray sources and compare the simulated dipole to that measured by Auger. We use the Fermi-LAT 4FGL-DR4 catalog as the basis for our simulations. We find that the resolved gamma-ray sources are insufficient to account for the population of sources producing UHECRs >8 EeV, primarily because of low density of nearby sources. We conclude that there must exist a population of UHECR sources that are gamma-ray dim.

Our assumption of a turbulent EGMF pervading the Universe is not realistic. Ideally, the distribution of magnetic fields in the large-scale structure of the Universe should be taken into account. In this case, cosmic voids would be the dominant component driving UHECR deflections (Alves Batista et al. 2017). By adopting a somewhat strong
EGMF ($B_{\text{RMS}}\sqrt{l_{\text{coh}}}$ = $10^{-13}$ G√Mpc), we are conservatively choosing a scenario that increases overall deflections and thus reduces the dipole amplitude, which strengthens our conclusions, in qualitative agreement with Dundović & Sigl (2019) and Hackstein et al. (2018).

As an additional test of the effect of EGMF strength on our simulated dipole, we also consider a larger EGMF strength of $B_{\text{RMS}}\sqrt{l_{\text{coh}}}$ = 0.16nG. While this value exceeds the upper bound on cosmological magnetic fields by Jedamzik & Saveliev (2019), it was obtained as the best-fit value in van Vliet et al. (2021), which uses simulations of UHECRs from nearby star-forming galaxies to reproduce the intermediate-scale anisotropy observed by Auger. In both energy flux and uniform weight scenarios, the larger magnetic field strength does not alter the dipole amplitude sufficiently to change the conclusion of our work. In the uniform weight case, the order of magnitude change in $B_{\text{RMS}}$ only lowers the dipole by a few percent. In the energy flux weight case, the larger EGMF strength actually raises the dipole by ~20%, this effect caused by the effect of EGMF on the balance between the brightest source, Markarian 421, and the other background sources.

While the JF12 model of the GMF has been found to fit many observable parameters well, the model is still not perfectly constrained, and it is expected that new GMF models will continue to improve our understanding of UHECR deflection and source searches in the future (Boulanger et al. 2018; Coleman et al. 2023). Other works (van Vliet et al. 2021; Bister & Farrar 2024) using the JF12 model make the assumption that variations in the GMF model can affect the amplitude and direction of the measured dipole somewhat, but that the overall trend of the dipole's behavior given a particular source distribution and EGMF strength would not change drastically. In the case of this work, while future models of the GMF could lower the evaluated dipole amplitude somewhat, a new GMF model would have to employ significantly larger magnetic field strengths to bring our evaluated dipole amplitude into agreement with that measured by Auger. The recent work Unger & Farrar (2024) provides updated models of the GMF, but no model in the new ensemble differs enough from JF12 to significantly lower the expected dipole amplitude.

In general, our finding is in agreement with previous works that have shown that it is necessary to have a very high density of nearby sources to produce a dipole amplitude as low as that observed by Auger (Ding et al. 2021; van Vliet et al. 2021; Bister & Farrar 2024; Abdul Halim et al. 2024). Given our assumptions about the EGMFs and GMFs, our result shows that the density of UHECR sources must be larger than that of the resolved gamma-ray sources.

Future neutrino telescopes (Ackermann et al. 2022) will seek to observe ultra-high-energy (UHE) neutrinos and to identify their sources. UHE neutrinos can be produced in UHECR sources (e.g., Fang et al. 2014), and their production is accompanied by that of UHE gamma rays. The UHE gamma rays may cascade down to GeV energies by the radiation background of the source or the extragalactic background lights. Our result implies that most UHECR accelerators might lack the conditions necessary for significant hadronic interactions to occur. Such sources may not only be gamma-ray dim, but also dim emitters of UHE neutrinos.

## Acknowledgments

The work of A.P. and K.F. is supported by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation. K.F. acknowledges support from National Science Foundation (PHY-2110821, PHY-2238916) and NASA (NMH211ZDA001N-Fermi). This work was supported by a grant from the Simons Foundation (00001470, KF). R.A.B. is funded by: the “la Caixa” Foundation (ID 100010434) and the European Union’s Horizon 2020 research and innovation program under the Marie Skodowska-Curie grant agreement No. 847648, fellowship code LCF/BQ/PI21/11830030; grants PID2021-125331NB-I00 and CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033, by “ERDF A way of making Europe,” and the MULTIDARK Project RED2022-134411-T.

Software: CRPRopa3.2 Alves Batista et al. (2022). We have made the code for converting CRPropa output to events observed by Auger publicly available at https://github.com/angelinapartenheimer/CRDipy.git.

## Appendix A
Spectrum at Earth and Flux Map in the Case of Gamma-Ray Flux Weighting

Here we present the flux map (Figure 5) and spectrum at Earth (Figure 6) of simulated cosmic rays in the case of gamma-ray flux weighting, which yields a simulated dipole amplitude of 74%. In this case, the UHECR dipole is dominated by Markarian 421, an especially bright and nearby source in the 4FGL-DR4 catalog. The effect of Markarian 421 is demonstrated by the bright spot in the center of the flux map. Because the UHECR signal is dominated by one source, the simulated UHECR spectrum deviates slightly from the Auger data at high energies. However, this small amount of disagreement cannot explain the order of magnitude in difference between the simulated dipole and the 6% dipole observed by Auger.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Flux map corresponding to the case of gamma-ray flux weighting. The flux map is obtained by taking the ratio of the smoothed count density and exposure map following Aab et al. (2017). The sky map is smoothed with a 45° top-hat function using the healpy.sphtfnc.beam2bl function.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Spectrum at Earth of simulated UHECRs from Fermi-LAT sources in the case of gamma-ray flux weighting. Because the signal is dominated by one source, there is some discrepancy between the simulated spectrum and Auger data at high energies; however, the disagreement in the spectra is not sufficient to explain the order of magnitude in difference between the simulated and observed dipole amplitude.</figcaption></figure>

## Appendix B
### Rayleigh Analysis

This section recapitulates the Rayleigh analysis method outlined in Aab et al. (2017).

The dipole is calculated in terms of its components perpendicular $d_⊥$ and parallel $d_z$ to the Earth's rotation. The complete dipole is given by

$$d = \sqrt{d_{\perp}^2 + d_z^2}. \qquad (B1)$$

To reconstruct the component of the dipole perpendicular to Earth's rotation $d_⊥$, we perform Rayleigh analysis in R.A. $α$. Applying the Galactic magnetic lens returns a density map of cosmic-ray events. Each position in the sky (R.A. $α_i$, decl. $δ_i$) has a corresponding weight $w_i$ corresponding to the event density at that position. The first-harmonic Fourier components are

$$a_{\alpha} = \frac{2}{N} \sum_{i=1}^{N_{pix}} w_i \cos \alpha_i, \quad b_{\alpha} = \frac{2}{N} \sum_{i=1}^{N_{pix}} w_i \sin \alpha_i, \qquad (B2)$$

where $N = \sum_{i=1}^{N_{pix}} w_i$ and $N_{pix}$ is the number of pixels in the sky map.

The amplitude and phase of the first harmonic of modulation are

$$r_{\alpha} = \sqrt{a_{\alpha}^2 + b_{\alpha}^2}, \quad \tan \phi_{\alpha} = \frac{b_{\alpha}}{a_{\alpha}}. \qquad (B3)$$

The dipole component $d_⊥$ is then given by

$$d_{\perp} \approx \frac{r_{\alpha}}{\langle \cos \delta \rangle}, \qquad (B4)$$

where $\langle \cos \delta \rangle = \frac{1}{N} \sum_{i=1}^{N_{pix}} \cos \delta_i$ is the mean cosine of decl. of detected events. The corresponding R.A. of the dipole is

$$\alpha_d = \phi_{\alpha}. \qquad (B5)$$

To reconstruct the dipole component parallel to the Earth's rotation $d_z$, we perform Rayleigh analysis in azimuthal angle $φ$.

This yields first-harmonic Fourier components $a_φ$ and $b_φ$. In this case, the dipole component $d_z$ is given by

$$d_z \approx \frac{b_{\phi}}{\cos l_{\text{obs}} \langle \sin \theta \rangle}, \qquad (B6)$$

where $l_{\text{obs}}$ is the latitude of the Pierre Auger observatory, and $\langle \sin \theta \rangle$ is the mean sign of event zenith angles. The decl. of the dipole is

$$\tan \delta_d = \frac{d_z}{d_{\perp}}. \qquad (B7)$$

## Appendix C
### Results from 3FHL and 4LAC-DR3 Sources

We test the dipole amplitude and direction of UHECRs coming from 3FHL or 4LAC-DR3 sources alone. The 3FHL catalog provides sources observed by Fermi-LAT in the range 10 GeV–2 TeV from 7 yr of data, and is the first Fermi-LAT catalog to use Pass 8 event-level analysis (Ajello et al. 2017). It contains 1556 sources. The 4LAC catalog gives data on the high-latitude active galactic nuclei on 50 MeV–1 TeV (Ajello et al. 2020), and it is a subset of the 4FGL-DR3 catalog. It contains 2863 sources. The 4LAC-DR3 catalog contains no Galactic sources; the 3FHL catalog contains some Galactic sources, and these are removed following the same method outlined in Section 2. This results in 197 sources being removed from 3FHL, 72 of which are unassociated. 3FHL has redshift data for 35% of its extragalactic sources; 4LAC-DR3 contains redshift data for 63% of its sources. The sources without redshift data provided are assigned random redshifts following the SFR distribution from Yüksel et al. (2008).

We perform separate simulations for UHECRs from the 3FHL and 4LAC-DR3 sources. In each case, both energy flux and uniform source injection weighting are tested. Results from the tests are shown in Table 2. In all cases, the dipole amplitude is significantly larger than observed by Auger.

<table><caption>**Table 2**
Dipole of UHECRs from Separate Simulations of 3FHL and 4LAC-DR3 Sources</caption><thead><tr><th>&lt;8 EeV</th><th>Dipole Amplitude (%)</th><th>α<sub>d</sub> (deg)</th><th>δ<sub>d</sub> (deg)</th></tr></thead><tbody><tr><td>3FHL flux weight</td><td>128</td><td>205</td><td>61</td></tr><tr><td>3FHL uniform weight</td><td>71</td><td>295</td><td>39</td></tr><tr><td>4LAC flux weight</td><td>77</td><td>253</td><td>70</td></tr><tr><td>4LAC uniform weight</td><td>31</td><td>303</td><td>-28</td></tr></tbody></table>

### ORCID iDs

Angelina Partenheimer https://orcid.org/0009-0007-2369-258X
Ke Fang https://orcid.org/0000-0002-5387-8138
Rafael Alves Batista https://orcid.org/0000-0003-2656-064X
Rogerio Menezes de Almeida https://orcid.org/0000-0003-3104-2724
References

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

Aab, A., Abreu, P., Aglietta, M., et al. 2018a, *ApJ*, 868, 4

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

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

Aab, A., Abreu, P., Aglietta, M., et al. 2020b, *PhRvD*, 102, 062005

Abbasi, R. U., Abe, M., Abu-Zayyad, T., et al. 2014, *ApJL*, 790, L21

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, *ApJS*, 247, 33

Abdul Halim, A., Abreu, P., Aglietta, M., et al. 2023, *PoS*, 444, 521

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

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

Ackermann, M., Bustamante, M., Lu, L., et al. 2022, *JHEAp*, 36, 55

Ajello, M., Angioni, R., Axelsson, M., et al. 2020, *ApJ*, 892, 105

Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, *ApJS*, 232, 18

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

Alves Batista, R., Dundovic, A., Erdmann, M., et al. 2016, *JCAP*, 2016, 038

Alves Batista, R., & Saveliev, A. 2020, *ApJL*, 902, L11

Alves Batista, R., Shin, M.-S., Devriendt, J., Semikoz, D., & Sigl, G. 2017, *PhRvD*, 96, 023010

Alves Batista, R., Tjus, J. B., Dörner, J., et al. 2022, *JCAP*, 2022, 035

Amerio, A., Cuoco, A., & Fornengo, N. 2023, *JCAP*, 2023, 029

Ballet, J., Bruel, P., Burnett, T. H., Lott, B. & The Fermi-LAT Collaboration 2023, arXiv:2307.12546

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

Boulanger, F., Enßlin, T., Fletcher, A., et al. 2018, *JCAP*, 2018, 049

Coleman, A., Eser, J., Mayotte, E., et al. 2023, *APh*, 147, 102794

di Matteo, A., & Tinyakov, P. 2018, *MNRAS*, 476, 715

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

Dundović, A., & Sigl, G. 2019, *JCAP*, 2019, 018

Fang, K., Kotera, K., Murase, K., & Olinto, A. V. 2014, *PhRvD*, 90, 103005

Fermi-LAT Collaboration 2018, *ApJS*, 237, 32

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

Hackstein, S., Vazza, F., Brüggen, M., Sorce, J. G., & Gottlöber, S. 2018, *MNRAS*, 475, 2519

Jansson, R., & Farrar, G. R. 2012, *ApJL*, 761, L11

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

Kotera, K., & Olinto, A. V. 2011, *ARA&A*, 49, 119

MAGIC Collaboration, Acciari, V. A., Agudo, I., et al. 2023, *A&A*, 670, A145

Unger, M., & Farrar, G. R. 2023, arXiv:2312.02645

van Vliet, A., Palladino, A., Taylor, A., & Winter, W. 2021, *MNRAS*, 510, 1289

Waxman, E., & Miralda-Escude, J. 1996, *ApJL*, 472, L89

Yüksel, H., Kistler, M. D., Beacom, J. F., & Hopkins, A. M. 2008, *ApJL*, 683, L5

Yushkov, A. 2019, *ICRC (Madison, WI)*, 36, 482