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

Ultra-high-energy cosmic rays, 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.


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(Aab et al. , 2018a(Aab et al. , 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 TA 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 CR-Propa3 (Abdollahi et al. 2020;Ajello et al. 2017Ajello et al. , 2020;;Alves Batista et al. 2016, 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.
The Large Area Telescope (LAT) aboard 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-year Source Catalog (4FGL-DR4; Ballet et al. 2023;Abdollahi et al. 2020), 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 7,195 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 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 highlatitude 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 skymap 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 de- scribed 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.

SIMULATION
Without loss of generality, we consider two source weighting schemes.In the first scenario, we assume that the 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. 2022Batista et al. , 2016)).We use an injection spectrum following a power law and exponential cutoff at maximum rigidity: We adopt a spectral index γ = −1.82 and maximum rigidity R max = 10 18.15 V. We adopt luminosity fractions I p = 0.620, I He = 0.233, I N = 0.533, I Si = 0.148, and I 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: root mean square value of the magnetic fields B RMS , and maximum correlation length l coh .In particular, we adopt 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 λ is given by (Waxman & Miralda-Escude 1996): The angular distribution of cosmic rays after averaging over time is given by a simple Gaussian: 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 ϕ 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, τ ≈ Dθ 2 s /(4c), corresponding to a distance D add = Dθ 2 s /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, 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 pc 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 ∈ [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 skymaps of UHECRs using CRPropa's ParticleMapsContainer feature, which yields a normalized skymap 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.

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.

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. 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 skymap is smoothed with a 45 • top-hat function using the healpy.sphtfnc.beam2blfunction.
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 ten 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 ten times larger than that measured by Auger.In this case, the dipole is dominated by 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 sig-nificantly 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 gammaray 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.

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 dN dS ∝ S −2 below the 4FGL-DR4 sensitivity of 2 × 10 −10 ph cm −2 s −1 .Following this distribution, we generate a population of unresolved sources extending to the minimum luminosity of S min = 5 × 10 −12 ph cm −2 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 dNγ dEγ ∝ E −bγ γ with b γ = 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 × 10 −4 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 2023) 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.

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 gammaray 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 extragalactic magnetic field (EGMF) pervading the universe is not realistic.Ideally, the distribution of magnetic fields in the largescale 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 RMS √ l 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 RMS √ l 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 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 2023) 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 which 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 (van Vliet et al. 2021;Ding et al. 2021;Bister & Farrar 2023;Abdul Halim et al. 2024).Given our assumptions about the extragalactic and Galactic magnetic fields, 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.We have made the code for converting CRPropa output to events observed by Auger publicly available at https: //github.com/angelinapartenheimer/CRDipy.gitBecause 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.
ysis in right ascension α.Applying the Galactic magnetic lens returns a density map of cosmic ray events.Each position in the sky (right ascension α i , declination δ i ) has a corresponding weight w i corresponding to the event density at that position.The first-harmonic Fourier components are: The dipole component d ⊥ is then given by: where ⟨cos δ⟩ = 1 N Npix i=1 cos δ i is the mean cosine of declinations of detected events.The corresponding right

Figure 2 .
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.

Figure 3 .
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).
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 followingAab et al. (2017).The skymap is smoothed with a 45 • top-hat function using the healpy.sphtfnc.beam2blfunction.

Figure 6 .
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.
N pix is the number of pixels in the skymap.The amplitude and phase of the first harmonic of modulation are: This choice is motivated by the upper bound on cosmological magnetic fields by Jedamzik & Saveliev (2019) (B RMS ≲ 10 −11 G).It is also in agreement with lower bounds derived from gamma-ray observations, which generally indicate B RMS ≳ 10 −17 G at coherence lengths l coh ≳ 100 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 non-existent or very weak (e.g., Alves Batista & Saveliev 2020).