# Constraints on UHECR sources and extragalactic magnetic fields from
directional anisotropies

Teresa Bister<EMAIL_ADDRESS>Institute for Mathematics, Astrophysics and
Particle Physics, Radboud University Nijmegen, Nijmegen, The Netherlands
Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF),
SciencePark, Amsterdam, The Netherlands  Glennys R. Farrar<EMAIL_ADDRESS>Center
for Cosmology and Particle Physics, New York University, New York, NY 10003,
USA

###### Abstract

A dipole anisotropy in ultra-high-energy cosmic ray (UHECR) arrival
directions, of extragalactic origin, is now firmly established at energies
$E>8\,\mathrm{EeV}$. Furthermore, the UHECR angular power spectrum shows no
power at smaller angular scales than the dipole, apart from hints of possible
individual hot or warm spots for energy thresholds $\gtrsim$ 40 EeV. Here, we
exploit the magnitude of the dipole and the limits on smaller-scale
anisotropies to place constraints on two quantities: the extragalactic
magnetic field (EGMF) and the number density of UHECR sources or the
volumetric event rate if UHECR sources are transient. We also vary the bias
between the extragalactic matter and the UHECR source densities, reflecting
whether UHECR sources are preferentially found in over- or under-dense
regions, and find that little or no bias is favored.

We follow Ding et al. (2021) in using the Cosmic Flows 2 density distribution
of the local universe as our baseline distribution of UHECR sources, but we
improve and extend that work by employing an accurate and self-consistent
treatment of interactions and energy losses during propagation. Deflections in
the Galactic magnetic field are treated using both the full JF12 magnetic
field model, with random as well as coherent components, or just the coherent
part, to bracket the impact of the GMF on the dipole anisotropy. This Large
Scale Structure (LSS) model gives good agreement with both the direction and
magnitude of the measured dipole anisotropy and forms the basis for
simulations of discrete sources and the inclusion of EGMF effects.

Cosmic ray astronomy (324); Cosmic ray sources (328); Cosmic rays (329);
Extragalactic magnetic fields (507); Milky Way magnetic fields (1057); Ultra-
high-energy cosmic radiation (1733); Cosmic anisotropy (316); Large-scale
structure of the universe (902)

††software: CRPropa3 (Alves Batista et al., 2016), healpy (Zonca et al.,
2019), numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), astrotools
(Bister et al.)

##  1 Introduction

Ultra-high-energy cosmic rays (UHECRs), characterized by energies exceeding
$10^{18}\,\mathrm{eV}$, continue to challenge our understanding of
astrophysical processes in the Universe. Large numbers of UHECRs have been
detected by the giant Pierre Auger Observatory (The Pierre Auger
Collaboration, 2015) (Auger) in the southern hemisphere and the Telescope
Array (The Telescope Array Collaboration, 2008) (TA) in the northern
hemisphere, leading to robust determinations of the UHECR energy spectrum and
arrival directions, and approximate information on the mass composition.

A major obstacle to finding the sources of UHECRs arises from their
interactions with magnetic fields during their propagation from sources to
Earth, which – given the uncertainty on the UHECR charge, $Z$, and the
intervening magnetic fields – impedes straightforward source identification.
Understanding UHECR trajectories demands a good knowledge of the Galactic
magnetic field (GMF) within our Milky Way. For that, several models have been
in the literature; see, e.g., the review by Jaffe (2019). A comprehensive new
suite of models encompassing the range of possible coherent GMF models
consistent with the latest data has very recently been released (Unger &
Farrar (2023); UF23 below). In this paper, we use the deflections predicted
with the classic JF12 model including random component (Jansson & Farrar,
2012a), and bracket the sensitivity to the random field by using the JF12
coherent field only. Unger & Farrar (2023) checked that deflections in the
JF12 coherent field generally fall within the range found for the UF23 suite
of updated models, for the rigidities tested, so the results of the present
study should be qualitatively robust.

In addition to the GMF, the extragalactic magnetic field (EGMF) plays a role
at some level in UHECR propagation. The strength, spatial distribution, and
even origin of the EGMF is unclear both theoretically and observationally.
Limits on parameters of the EGMF exist from theoretical expectations and
observations of the CMB or $\gamma$-rays (see review Durrer & Neronov (2013)),
but large parts of the parameter space are still unconstrained. One of the
accomplishments of the present work is placing new constraints on key
parameters of the EGMF through a combination of its field strength and
coherence length.

The observation of large-scale anisotropies in the arrival directions of
UHECRs allows for valuable insights into the distributions of sources, and at
the same time into the interactions with cosmic magnetic fields. The Pierre
Auger collaboration made the important discovery of a large-scale dipole
anisotropy of the UHECR arrival directions above 8 EeV (The Pierre Auger
Collaboration, 2017a, 2018), with a magnitude of $7.3\%$ above 8 EeV and a
significance of $6.9\sigma$ (Golup, G. for the Pierre Auger Collaboration,
2023). The Telescope Array data supports these findings, but without being
able to significantly differentiate from isotropy due to limited statistics
(The Telescope Array Collaboration, 2020). The dipole strength increases with
the energy threshold. Its direction is well away from the Galactic center,
adding to evidence for an extragalactic origin of UHECRs at these energies.
Another important Auger result is that all higher multipole moments are
compatible with isotropy (de Almeida, R. for the Pierre Auger Collaboration,
2021). This is also corroborated by results of the combined working group
between Auger and TA (Caccianiga, L. for the Pierre Auger and Telescope Array
Collaborations, 2023).

In the absence of an obvious dominant source, a natural hypothesis is that
UHECR sources follow the local matter density, and thereby the large-scale
structure (LSS) of the Universe. This scenario was investigated in e.g. Harari
et al. (2010, 2015); Globus & Piran (2017); Globus et al. (2019); The Pierre
Auger Collaboration (2017a); Ding et al. (2021) with promising results
indicating that the dipole amplitude and direction may be due to the LSS.
Here, we make significant technical improvements on the work of Ding et al.
(2021). Our results corroborate the conclusion that the LSS naturally can
explain the observed dipole. We then exploit the dipole and higher multipole
anisotropy observations to draw new and important conclusions about UHECR
sources and extragalactic magnetic fields. Specifically, we place limits on
the _bias_ between UHECR sources and the LSS (whether UHECR sources populate
predominantly matter-dense or underdense regions), and we constrain
extragalactic magnetic fields. This enables us to place improved constraints
on the source density, for steady sources, or the volumetric rate if sources
are transients. These constraints allow the possible classes of UHECR sources
to be further delimited.

Some of these ideas were investigated in Allard et al. (2022), but without a
fit of the injection parameters and with a different treatment of propagation
effects. A distinction between Allard et al. (2022) and our analysis is the
different approaches to modeling the UHECR source distribution. Determination
of the local galaxy or matter distribution is ultimately strongly reliant on
the 2MRS flux-limited catalog, which excludes 8∘ from the Galactic plane (5∘
in the outer Galaxy). To avoid distance-dependent selection biases, users must
create a volume-limited sub-catalog. To do that to $176\,\mathrm{Mpc}$, given
the 2MRS flux limit, requires restricting to such bright sources that the
number density is $1.1\times 10^{-4}\,{\rm Mpc}^{-3}$ (Allard et al., 2022),
just 1% the density of ordinary galaxies (Conselice et al., 2016). This biases
the source distribution to very bright galaxies and limits the ability to make
independent samples of sources except for low source densities, which they
address by some procedure to generate synthetic galaxy catalogs at higher
density. This also requires that large peculiar velocities, leading to noise
in the determination of galaxy distances, be corrected which is challenging.
Instead of trying to do those things ourselves, we adopt the CosmicFlows2
(Hoffman et al., 2018) 3D matter density, which was determined by self-
consistently modeling galaxy peculiar velocities as well as distances, angular
positions, and redshifts, including very careful analysis of all the available
distance measures for each galaxy used in their modeling. To our mind, this is
the optimal conceptual approach to finding the matter density distribution in
the angular region observed by 2MRS, and moreover thanks to self-consistent
modeling of the potential in which observed galaxies move, and inclusion of
more observations, the mass distribution in the region behind the Galactic
plane is also approximately described. We find that approximating the UHECR
source density by the CosmicFlows matter density without a bias gives a good
fit to the observations, and that restricting to regions of high mass density
leads to significantly worse agreement with the dipole.

Our paper is structured as follows. First, the analysis pipeline is introduced
in Sec. 2, including the modeling of the source setup, propagation, magnetic
fields, and the fit procedure. The results for the baseline case are presented
in Sec. 3, and the influence of possible model uncertainty is elaborated in
Appendix A. In Sec. 4, we investigate the possibility of a bias between the
UHECR source distribution and the matter density distribution. Then, the
constraining power of the model regarding the extragalactic magnetic fields as
well as the source number density is investigated in Sec. 5, and astrophysical
implications of the constraints are presented. Finally, we conclude with a
summary in Sec. 6.

##  2 Analysis pipeline: modeling UHECRs from sources to Earth

In this section, the analysis pipeline is described, i.e. the injection,
propagation, and the deflections in magnetic fields. The fit method and
likelihood function are presented in Sec. 2.3.

###  2.1 Source distribution and injection

The baseline distribution of UHECR sources is taken to follow the
(predominantly dark) matter density of the nearby universe, as determined in
by CosmicFlows2 (Hoffman et al., 2018) which is based on self-consistently
combining peculiar galaxy velocities with redshifts and angular positions. We
use the CosmicFlows densities between 0 and 360 Mpc in shells of 1.3 Mpc with
an angular binning using healpy (Zonca et al., 2019) with nside=32 (leading to
roughly 12000 voxels per shell). At larger distances, we assume a homogeneous
distribution of sources up to 5000 Mpc using a linear extrapolation of the
density. We take the injection of UHECRs to be common to all sources and model
the injection rate (number of nuclei injected per unit of energy, volume, and
time) to follow a Peters cycle:

$Q(E,A)\propto\ a_{A}\ E^{\gamma}\ f_{\mathrm{cut}}(E,Z_{A}R_{\mathrm{cut}}).$
(1)

It is characterized by a spectral index $\gamma$ and a cutoff at the maximum
rigidity $R_{\mathrm{cut}}$, described by a variable cutoff function
$f_{\mathrm{cut}}(x=\frac{E}{Z_{A}R_{\mathrm{cut}}})$. The cutoff function is
either a broken exponential (b.e.) of the form
$f_{\mathrm{cut}}^{\mathrm{b.e.}}(x)=\exp(1-x)$ (that sets in for $x>1$) to
compare to Auger results (The Pierre Auger Collaboration, 2017b), or a
continuous hyperbolic secans function with variable softness
$f_{\mathrm{cut}}^{\mathrm{sec}}(x,k)=1/\text{cosh}(x^{k})$ with $k=1,2,3$ as
in Mollerach & Roulet (2019). We inject and track five representative elements
with charge numbers $Z_{A}$ and contributions below the cutoff $a_{A}$: H, He,
N, Si, and Fe. Instead of the fraction $a_{A}$ defined below the cutoff, we
quote below the integral fractions $I_{A}=a_{A}\
Z_{A}^{2-\gamma}/(\sum_{A}a_{A}\ Z_{A}^{2-\gamma})$ as these represent the
total emission better.

###  2.2 Propagation and cosmic magnetic fields

For the propagation, we use a database of 1-dimensional CRPropa3 (Alves
Batista et al., 2016) simulations for which all available interactions are
considered. For the extragalactic background light, the Gilmore model (Gilmore
et al., 2012) is used, and for the photodisintegration the TALYS model (Koning
et al., 2005). The interactions with background photon fields lead to a
shrinking propagation horizon for increasing energies, as can be seen in Fig.
1.

Figure 1: Map in Galactic coordinates showing the mean distance from which a
UHECR originates in the baseline model, for various energy bins or energy
thresholds. Both the anisotropic extragalactic matter distribution as well as
the effects of propagation are visible in the maps. With increasing energy,
the overall propagation horizon shrinks. At the highest energies, many CRs
come from nearby, predominantly from the Virgo cluster (compare to Fig. 2).

As in Ding et al. (2021), we calculate illumination maps $I$. These encode the
UHECR flux at the edge of our Galaxy induced by the local anisotropic source
distribution and more distant homogeneous distribution, shaped by propagation
effects. The illumination maps for our best-fit model, assuming zero EGMF, are
displayed in Fig. 2, for different distance bins as well as one cumulative
bin. Evidently, the illumination map is peaked in the Galactic North, where
the Virgo cluster and several other massive clusters reside.

Figure 2: Logarithmic illumination maps (see text) for different distance
shells, and a cumulative one (lower right), for the best-fit model and $E>8$
EeV. The direction of the full-sky dipole component of the map is indicated
with a grey marker. For the closest distance bin, $0-40$ Mpc, the influence of
the Virgo cluster in the Galactic north is clearly visible. For the further
away distances, mainly Perseus-Pisces and the Great Attractor become important
at $40-80$ Mpc, and Coma at $80-120$ Mpc.

The local arrival direction distribution can be calculated from the
illumination maps by using a lens of the Galactic magnetic field (GMF). As a
baseline case, we use the same high-resolution backtracking simulations as
Ding et al. (2021) of the JF12 model (Jansson & Farrar, 2012a, b) with a
coherence length of 30 pc for the random component. As pointed out in Farrar
(2014), the turbulent component of the GMF is likely overestimated in the
original JF12 random model (Jansson & Farrar, 2012b) since JF12 was based on
the original WMAP synchrotron intensity map which has subsequently been
revised downward due to better modeling of dust contributions (The Planck
collaboration, 2016). To bracket the effect of the random field, we therefore
also consider the case with only the JF12 regular magnetic field without a
random component, denoted as JF12 reg.

Following the results from Ding et al. (2021) who find that very weak
extragalactic magnetic field (EGMF) strengths lead to the best agreement with
the data, we neglect the extragalactic magnetic field in the baseline case.
This is backed up by evidence for small field strengths in voids (The Planck
collaboration, 2016; Hackstein et al., 2016). In Sec. 5.1 we constrain the
maximum allowed strength of a turbulent EGMF. The effect of magnetic
confinement by large magnetic field strengths in denser galaxy clusters is
investigated in the context of bias in Sec. 4.

###  2.3 Fit of model parameters

After propagation and magnetic field deflections, the arriving (detected) CRs
are binned into a matrix whose elements are
$(E_{\mathrm{det}},A_{\mathrm{det}},D,\mathrm{pixel})$. Here, $\mathrm{pixel}$
labels the healpix sky-direction and $D$ the propagation distance, binned as
described above. $E_{\mathrm{det}}$ labels bins in energy of width
$10^{0.1}\,\mathrm{eV}$, between $10^{18.0}\,\mathrm{eV}=1\,\mathrm{EeV}$ and
$10^{20.4}\,\mathrm{eV}\approx 250\,\mathrm{EeV}$. The masses
$A_{\mathrm{det}}$ are binned as $A_{\mathrm{det}}\in(1,\ 2\text{--}4,\
5\text{--}22,\ 23\text{--}38,\ $>39), to compare to the fitted Auger
composition (The Pierre Auger Collaboration, 2023a, b).

We use a Poissonian likelihood $\mathcal{L}_{E}$ as given in The Pierre Auger
Collaboration (2017b) to compare the modeled energy spectrum above
$10^{18.9}\,\mathrm{EeV}\approx 8\,\mathrm{EeV}$ to the measured one. (This is
the same energy threshold as for the dipole components, see below). We use the
unfolded spectrum from (The Pierre Auger Collaboration, 2020), so that no
detector effects on the energy spectrum have to be considered in the model.
Below $10^{18.9}\,\mathrm{EeV}$, we ensure that the model does not overshoot
the data by a one-sided Poissonian likelihood.

We use the peak of the shower profile, $X_{\mathrm{max}}$, as the mass
estimator. EPOS-LHC (Pierog et al., 2013) is used as the baseline hadronic
interaction model and Sibyll2.3d (Riehn et al., 2020) as an alternative. The
effects induced by the detector during measurement (The Pierre Auger
Collaboration, 2014) are included in the modeling of $X_{\mathrm{max}}$. As in
(The Pierre Auger Collaboration, 2023b), we consider a possible shift of the
$X_{\mathrm{max}}$ scale as a nuisance parameter in the model. That parameter
$\nu_{X_{\mathrm{max}}}$ is given in units of standard scores that refer to
multiples of the experimental systematic uncertainty, so that the likelihood
function $\mathcal{L}_{\mathrm{syst}}$ is a Gaussian (see The Pierre Auger
Collaboration (2023b) for details). Finally, the modeled shower depth
distributions binned in energy $>10^{18.9}\,\mathrm{EeV}$ are compared to the
data from Yushkov, A. for the Pierre Auger Collaboration (2019) via a
Multinomial likelihood $\mathcal{L}_{X_{\mathrm{max}}}$ as in The Pierre Auger
Collaboration (2017b).

To compare to the measurements of the large-scale anisotropy, the directional
exposure of the Observatory is taken into account. Afterwards, we acquire the
modeled arrival directions in energy bins
$E_{\mathrm{bins}}=(8-16\,\mathrm{EeV},16-32\,\mathrm{EeV},>32\,\mathrm{EeV})$,
for which we calculate the 3D dipole components $d_{x}$, $d_{y}$, and $d_{z}$
considering the effect of the limited exposure ((Deligny & Salamida, 2013);
see also K-inverse method explained in Ding et al. (2019)). The dipole
components $d_{e,i\in(x,y,z)}$ for each of the three energy bins are then
compared to the measured $d_{e,i\in(x,y,z)}^{\mathrm{data}}$ (including the
uncertainties $\sigma_{e,i\in(x,y,z)}^{\mathrm{data}}$) from Auger (de
Almeida, R. for the Pierre Auger Collaboration, 2021) using a Gaussian
likelihood function:

$\mathcal{L}_{d}=\prod_{e\in
E_{\mathrm{bins}}}\prod_{i\in(x,y,z)}\frac{1}{\sqrt{2\pi}\sigma_{e,i}^{\mathrm{data}}}e^{-\frac{1}{2}\Big{(}\frac{d_{e,i}-d_{e,i}^{\mathrm{data}}}{\sigma_{e,i}^{\mathrm{data}}}\Big{)}^{2}}.$
(2)

That way, the energy evolution of both the direction and amplitude of the
dipole are taken into account.

Finally, all parts of the log-likelihood are added to get the total one:

$\log\mathcal{L}_{\mathrm{tot}}=\log\mathcal{L}_{E}+\log\mathcal{L}_{X_{\mathrm{max}}}+\log\mathcal{L}_{d}+\log\mathcal{L}_{\mathrm{syst}}$
(3)

##  3 Fit results

###  3.1 Composition and spectrum

Figure 3: Energy spectrum (upper) and first two moments of the
$X_{\mathrm{max}}$ distributions (lower) for the baseline model. The grey area
below 8 EeV is not included in the fit. The measured $X_{\mathrm{max}}$
distributions denoted by black markers have been moved by the best-fit
$X_{\mathrm{max}}$ scale shift (see text).

Figure 4: Frequency (a.u.) of arriving rigidities on Earth for the baseline
model for two energy thresholds, $E>8$ EV (left) and $E>32$ EV (right). The
multi-peak structure is due to our coarse binning in mass following Auger.

In the baseline case (CosmicFlows2 as source proxy, EPOS-LHC as hadronic
interaction model, JF12-reg+rand as GMF model, broken exponential at
injection), the model describes the measured spectrum and $X_{\mathrm{max}}$
data at Earth well, as depicted in Fig. 3. The corresponding best-fit
parameters and values of the likelihood are provided in Table 2 in the
Appendix. The best-fit injection parameters are close to recent results
(obtained with a different energy threshold) by Auger, using a homogeneous
source model (The Pierre Auger Collaboration, 2023b, 2017b). The best-fit
spectrum is quite hard, with spectral index $\gamma=-1.2$, the maximum
rigidity is $R_{\mathrm{cut}}=10^{18.2}$ V, and the composition is nitrogen-
dominated. For the $X_{\mathrm{max}}$ scale, a best-fit shift of
$\nu_{X_{\mathrm{max}}}=-0.9\sigma$ is found. This shift of the data (as
displayed in Fig. 3) leads to a better likelihood value, and thus a better
agreement of model and data, than no shift. It translates to a heavier
composition on Earth than without a shift, in agreement with the analysis by
Auger (The Pierre Auger Collaboration, 2023b).

Because the best-fit injected spectrum is very hard, the overlap of different
mass groups is small, as visible in Fig. 3. Consequently, heavier and heavier
mass groups dominate with rising energy. This parallel increase of the charge
$Z$ with the energy $E$ leads to a relatively constant rigidity $R=E/Z$ that
does not evolve with the energy. The rigidity spectra above two energy
thresholds, $E>8$ EeV and $E>32$ EeV, are depicted in Fig. 4. The mean
rigidities are 3.5 EV and 4.5 EV for the respective energy thresholds, and
almost all arriving UHECRs have rigidities between 1 EV and 5 EV111For
reference, $\log_{10}(5\times 10^{18})=18.7$. For the model using Sibyll
instead of EPOS-LHC as hadronic interaction model, (see sec A of the Appendix)
the mean rigidities for the two energy thresholds are 3.0 EV and 4.1 EV.. Note
that the multi-peak structure in the rigidity spectrum is because of the
coarse binning in mass following Auger. We tested that the effect of this on
the analysis is negligible by comparing to a model using the smoother mass
distribution from Muzio et al. (2019).

|  Illumination
---
|  flux, baseline, no TH
---
|  flux, baseline, TH=$45^{\circ}$
---
|  flux, JF12 reg, no TH
---
|  flux, JF12 reg, TH=$45^{\circ}$
---
Figure 5: Baseline model illumination maps $I$ (upper row) and arrival maps
$J$, without (second row) and with (third row) $45^{\circ}$-tophat (TH)
smearing. The columns show different energy thresholds or bins; see text below
each figure. The last two rows are for the JF12 regular-only GMF model with no
random component. All maps show the ratio of flux to the map mean, in Galactic
coordinates. The black point and circles indicate the measured dipole with
$68\%$ and $90\%$ domains from Auger (de Almeida, R. for the Pierre Auger
Collaboration, 2021). The pink star and circle show the direction and $68\%$
uncertainty of the model dipole, including the Observatory exposure; the grey
marker indicates the full-sky dipole direction in the LSS model.

###  3.2 Dipole predictions

Due to the energy-independence of the rigidity, as seen in Fig. 4, it is not
surprising that the observed dipole direction stays relatively constant for
different energy bins, both for observations and predictions, as can be seen
in the last four rows of Fig. 5. The Auger dipole is shown in black and the
model predictions, including the effect of exposure, in pink. At the highest
energies, the horizon is relatively close, where sources toward the Galactic
North Pole dominate the flux. This moves the dipole direction, even though the
rigidity is not much different. This horizon effect can be seen comparing the
illumination maps shown in the upper row of Fig. 5 (see also Fig. 2).

The model prediction of the dipole direction is depicted including its 68%
C.L. statistical uncertainty in Fig. 5. To estimate this statistical
uncertainty, we created $5,000$ realizations of $N_{\mathrm{events}}=44,000$
events (as in the Auger data set (de Almeida, R. for the Pierre Auger
Collaboration, 2021)), following the model. Earlier studies (Ding et al.,
2021; Allard et al., 2022) did not provide an estimation of the statistical
uncertainties in their dipole direction and amplitude predictions.

Note that based on the LSS model, no strong excesses are expected in the
Telescope Array field of view, indicating that the observed hotspot (Kim, J.
for the Telescope Array collaboration, 2023), if statistically significant,
should arise from a local source rather than the LSS. Note also that the
speculated correlation of an anisotropy with the Perseus-Pisces supercluster
(Kim, J. for the Telescope Array collaboration, 2023) is not observed in the
UHECR arrival maps in Fig. 5: even though there is an overdensity in the
illumination map from the supercluster (Fig. 2), it is displaced or dissolved
completely by the GMF.

Figure 6 shows the dipole amplitude. One sees that it increases with the
energy. This is due to a combination of the shrinking propagation horizon
(inducing greater anisotropy in the illumination maps) and to some extent a
decreasing diffusion in the GMF for higher-rigidity particles.

Both the predicted amplitude and direction of the dipole are very close to the
measurements for the highest energy bin >32 EeV. For smaller energies, the
dipole amplitude is well-bracketed by the baseline model using the JF12
coherent plus random field and the regular-field only model, shown with filled
and open stars respectively. However, for both models, the direction is
somewhat outside the 2$\sigma$ region. An exciting prospect is that the dipole
direction can be used to identify a preferred subset within the new suite of
models provided by Unger & Farrar (2023), designed to encompass the present
uncertainties in the GMF.

Figure 6: Dipole amplitude comparing to the Auger data (black markers, called
$d$ in Table 1 of (de Almeida, R. for the Pierre Auger Collaboration, 2021)).
The grey squares are the full-sky dipole moments of the illumination map (from
Fig. 5, upper row) and grey diamonds are the full-sky model prediction. The
pink stars are the model predictions including the effect of exposure, with
error bars showing the 68% C.L. statistical uncertainty due to the limited
number of events as discussed in the text, filled for the JF12 reg+rand
(baseline) model and unfilled for JF12-reg. The data points are positioned at
the energy-center of the data set for each energy threshold, with model
markers offset on the x-axis for better visibility.

###  3.3 Model variations

We tested the influence of the dipole (compared to spectrum and composition)
on the overall fit by removing the dipole likelihood part from Eq. 3. In that
case, where only energy spectrum and composition are considered, the fit
parameters as well as the likelihood values are quite similar to the baseline
case. The values are given in Table 2 in the Appendix. Clearly, the fit is
predominantly constrained by the spectrum and composition, and is not
distorted in order to fit the dipole components. This is confirmed by the fact
that the fit parameters do not differ much between the JF12-reg+rand
(baseline) and the JF12-reg cases, as also seen in Table 2.

The effects of the most important model uncertainties - the shape of the
injection cutoff and the uncertainty in the composition due to the hadronic
air shower modeling - are investigated in Sec. A of the Appendix, with the
respective fit parameters and likelihood values given in Table 2. It is
evident that the predictions of the model are independent of such systematic
uncertainties.

Thus, in the following, we use the baseline LSS model to obtain further
insights into the influence of cosmic magnetic fields, and also the number
density and bias of sources.

##  4 Production of UHECRs in high- and low-matter-density regions

In this section, we investigate the effect of a possible bias in the relation
between UHECR sources and the underlying LSS distribution. Such a bias could
arise, for example, if UHECR sources are predominantly found only in matter-
dense regions, as is the case for example for the most energetic active
galactic nuclei (AGNs) (Hale et al., 2017), or if UHECRs are accelerated in
the large scale shocks of massive galaxy clusters (Blandford et al., 2018);
see also Oikonomou et al. (2015) for exemplary simulations with varying bias.
On the other hand, a scenario where UHECRs are accelerated everywhere, but the
stronger magnetic fields inside matter-dense regions inhibit their escape
(e.g., Condorelli et al. (2023)), would lead to a bias towards low-density
regions.

|  baseline
---
|  cut $\rho<0.5\overline{\rho}$
---
|  cut $\rho<1.5\overline{\rho}$
---
|  cut $\rho>4\overline{\rho}$
---
Figure 7: Example illumination (left) and arrival (right) maps illustrating
the effect of bias. First row, the baseline case with no bias; second row, a
small bias away from low-density regions with
{$\rho_{\mathrm{min}}/\overline{\rho}=0.5$,
$\rho_{\mathrm{max}}/\overline{\rho}=\infty$}; third row, a large bias away
from low-density regions with {$\rho_{\mathrm{min}}/\overline{\rho}=1.5$,
$\rho_{\mathrm{max}}/\overline{\rho}=\infty$} (note different color bar
normalization!); and last row a small bias away from high density regions
{$\rho_{\mathrm{min}}/\overline{\rho}=0$,
$\rho_{\mathrm{max}}/\overline{\rho}=4.0$}. In comparison to the bias-free
illumination map (upper row), it is visible how the cut on high-density
regions partly removes the nearby Virgo cluster in the North. The cuts
excluding the low-density regions lead to more anisotropy in the illumination
and thereby to a larger dipole amplitude. For the extreme cut favoring the
high-density regions, the dipole direction is displaced more towards the
Galactic North than in the baseline case.

We test for both possible directions of bias by excluding regions of the
source distribution in the fit, if they either exceed a maximum density
$\rho_{\mathrm{max}}$, or are lower than a minimum density
$\rho_{\mathrm{min}}$ (where the mean density is denoted as
$\overline{\rho}$). Example illumination maps and corresponding flux maps are
depicted in Fig. 7, visualizing the impact of bias. Note that a bias towards
high-density regions leads to a shift of the predicted dipole direction
towards the Galactic North. This could explain slight differences in the
predicted dipole direction between our work and Allard et al. (2022), who use
a volume-limited, inevitably biased 2MRS-based catalog (c.f., discussion in
the Introduction).

The result of a scan over bias using multiple combinations of
$\rho_{\mathrm{min}}$ and $\rho_{\mathrm{max}}$ is shown in Fig. 8. It can be
seen that the likelihood decreases for almost any value of the cuts. Without
the high-density regions, the sky becomes more isotropic. Then, the dipole
direction is not represented well anymore, especially as even the conservative
cut $\rho_{\mathrm{max}}=4\overline{\rho}$ already removes parts of the nearby
Virgo cluster in the Galactic north (see Fig. 7, left), which is important to
reproduce the observed dipole. Removing low- and intermediate-density regions
with $\rho\lesssim\overline{\rho}$ leads to too large anisotropies,
overshooting the dipole amplitude. Only a cut removing the very low-density
regions below $0.5\,\overline{\rho}$ leads to a slight, non-significant
improvement of the likelihood. This study suggests that UHECR sources reside
in both high- and medium-density regions, with the CosmicFlows matter density
being a good proxy for the UHECR source distribution. No definite conclusion
can at this moment be drawn regarding low-density regions.

Figure 8: Influence of a bias between UHECR source distribution and LSS.
Regions with densities below $\rho_{\mathrm{min}}$ or above
$\rho_{\mathrm{max}}$ are excluded from the UHECR injection during the fit.
Colors and numbers indicate the likelihood ratio
($\log\mathcal{L}_{\mathrm{tot}}-\log\mathcal{L}_{\mathrm{tot}}(x)$) compared
to the baseline model without cuts (marked with white x).

##  5 Constraining the EGMF and the source number density

The observation of the dipole offers a unique opportunity to constrain the
intertwined relationship between UHECR sources and cosmic magnetic fields. The
dipole amplitude is heavily influenced by several effects, most importantly
the anisotropy in the source distribution and also the damping due to
turbulent magnetic fields. Within the LSS model, the source anisotropy is most
significantly influenced by the source number density and to some extent the
bias discussed in the previous section.

In this section, we investigate the effects of the EGMF and source number
density on the dipole. We begin by introducing the effect of the extragalactic
magnetic field (EGMF) on the dipole amplitude, for the continuous LSS source
model. Then, we explore the influence of different possible number densities
of UHECR sources on the dipole and higher multipole moments for the case of no
EGMF. Finally, we study a combined scenario with both EGMF and varying source
density.

###  5.1 Constraining the extragalactic magnetic field

We have demonstrated above that the LSS source model with the baseline or
JF12-reg models of the Galactic magnetic field, brackets the observed dipole
amplitude in the absence of significant extragalactic magnetic field smearing.
As described in the Introduction, the EGMF is poorly constrained by
observations. It can also have a significant effect on the dipole amplitude,
as well as higher order multipoles, through a smoothing of the flux
distribution before entering the Milky Way. In this subsection, we derive EGMF
constraints for the case of continuous source density.

To investigate the influence of the EGMF on the UHECR arrival directions, we
consider the most important effect expected for UHECRs propagating in a
turbulent field in the non-resonant scattering regime: a Gaussian beam
widening222We do not consider additional flux suppression due to longer path
lengths of low-rigidity particles diffusing in the EGMF, because given the
maximum possible values of $B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$ (see
Fig. 11) the increased path length is negligible, from Eq. 6. A more detailed
discussion can be found in Harari et al. (2015); Mollerach & Roulet (2019)..
When a UHECR passes through a turbulent magnetic field and its Larmor radius,
$r_{\text{g}}$, is large compared to the coherence length $L_{c}$ of the
field, its direction of propagation evolves diffusively in a calculable way.
The general relationship for the RMS deflection is clear:
$\delta\theta\sim\sqrt{DL_{c}}/r_{\text{g}}$, but a review of the literature
shows a variety of inconsistent expressions for the coefficient; this is
partly due to imprecision in defining the spectrum of turbulence and what is
meant by $L_{c}$, and also because some authors do not concern themselves with
factors of $\mathcal{O}(1)$. We find the analysis of Achterberg et al. (1999)
to be particularly useful and explicit. We adopt a Kolmogorov spectrum of
turbulence, for which the diffusion coefficient is derived at the end of App.
A (Achterberg et al., 1999): $\mathcal{D}_{0}=0.227L_{c}/r_{\text{g}}^{2}$
where $L_{c}$ is $\ell_{c}$ in the notation of Achterberg et al. (1999). The
RMS scattering angle is then

$\delta\theta=2.9^{\circ}\frac{B}{\mathrm{nG}}\frac{10\,\mathrm{EV}}{E/Z}\frac{\sqrt{D\
L_{c}}}{\mathrm{Mpc}}=2.9^{\circ}\beta_{\rm
EGMF}\frac{10\,\mathrm{EV}}{E/Z}\sqrt{\frac{\overline{D}}{\mathrm{Mpc}}},$ (4)

where we have introduced in the second equality the combination $\beta_{\rm
EGMF}\equiv B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$ to isolate the quantity
our analysis can constrain.

The mean distance $\overline{D}$ is calculated for each energy bin of width
$\log_{10}(E/\mathrm{eV})=0.1$ (compare to Fig. 1). To simulate the influence
of the EGMF, we convolve the illumination map for each rigidity $R=E/Z$ with a
Gaussian of width $\delta\theta$. It is not necessary to do a refit of the
source injection parameters for each value of $\beta_{\mathrm{EGMF}}$ because,
as demonstrated above, the fit is governed by the global energy spectrum and
composition observables. This is also supported by the finding that the fit
parameters do not really differ for the JF12 reg and JF12 reg+rand (baseline)
cases as discussed above.

Figure 9: Influence of a turbulent EGMF on the dipole amplitude for different
values of $\beta_{\rm EGMF}\equiv B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$,
compared to the JF12 reg model without turbulent GMF part.

The dipole amplitudes for different values of $\beta_{\mathrm{EGMF}}$ are
shown in Fig. 9 for the JF12-reg model, i.e., no turbulent GMF. Using JF12-reg
rather than the baseline model minimizes spreading due to the turbulent GMF
and hence enables us to place an upper bound on the spreading due to the EGMF
in order not to reduce the magnitude of the dipole too much. It can be seen in
Fig. 9 that for values $\beta_{\mathrm{EGMF}}\gtrsim 0.2$ the dipole amplitude
becomes significantly smaller than seen in the Auger data at lower energies.
Thus, in the limit of continuous sources, the extragalactic magnetic field
should respect $\beta_{\mathrm{EGMF}}\leq 0.2$. Hence for a coherence length
of, say, $L_{c}=1\,\mathrm{Mpc}$, the RMS field strength should not be
stronger than $B=0.2$ nG .

###  5.2 Constraining the source number density

In this subsection we place constraints on the source number density, $n$,
setting the extragalactic magnetic field to zero. (The lower limit quoted from
blazars (Aharonian et al., 2023) is so much lower than relevant for UHECR
smearing we can take it to be zero for this purpose.) We use both the baseline
LSS model and the JF12-reg model, which bracket the amount of random
deflections due to the GMF. A value of $n$ which is incompatible with the
dipole magnitude and power spectrum for either case can be excluded. When the
random component of the GMF is better constrained in the future, the exclusion
region will improve. In the next section, Sec. 5.3, the more general problem
of finding the allowed region in the source density-EGMF plane will be
treated.

Reducing the source density from the essentially continuum CosmicFlows model
to a finite number density, $n$, produces an additional source of
stochasticity in addition to the number of events in the UHECR data set
already discussed. This shows up in both the amplitude and direction of the
dipole and, importantly, in the rest of the angular power spectrum. Auger (de
Almeida, R. for the Pierre Auger Collaboration, 2021), as well as the combined
Auger+TA working group (Caccianiga, L. for the Pierre Auger and Telescope
Array Collaborations, 2023), has found that apart from the dipole, all higher
multipoles are consistent with being in the 99% range of variation found in
isotropic data sets. This proves to be highly constraining.

A virtue of using the power spectrum rather than the dipole direction to
constrain the number density is the insensitivity of the power spectrum to the
exact coherent GMF configuration (Tinyakov & Urban, 2015). This means that we
obtain robust constraints in spite of the non-perfect reproduction of the
dipole direction at lower energies. How the source number density affects the
dipole direction, as well as the energy dependence of the dipole, is reported
in Appendix B.

To model the effect of different source densities, we draw – for each $n$
studied – 1000 explicit catalogs of sources from the continuous source
distribution. Concretely, for speed reasons, for each catalog the closest
$10,000$ sources are drawn and the continuous distribution is used from the
distance of the most distant explicit source outwards. We verified that this
procedure does not affect the analysis results. From each source catalog, we
draw $N_{\mathrm{event}}=44,000$ events from those sources, that being the
event statistic in the Auger analysis (de Almeida, R. for the Pierre Auger
Collaboration, 2021), as we did for the continuum model analysis in Sec. 3.
Additionally, we draw 6000 events following the TA exposure and statistics as
in the combined data set used by the Auger+TA working group (Caccianiga, L.
for the Pierre Auger and Telescope Array Collaborations, 2023) in order to
have full-sky coverage to calculate the power spectrum.

For each realization, we ask:
1) _Are all the higher multipoles within the 99% isotropic expectation, for
all $C_{l>1}$ and energy ranges?_
2) _Is the dipole amplitude large enough?_

We use the dipole from Auger, rather than the combined Auger-TA dipole,
because even though having full-sky coverage in principle would give a more
accurate result, in practice the systematic uncertainty in the relative
calibration between observatories means that the Auger-only dipole is more
precise Caccianiga, L. for the Pierre Auger and Telescope Array Collaborations
(2023), as described in sec. 2.3. A too-large dipole amplitude for a given
source density can be mitigated by EGMF smoothing, so we only require that the
dipole amplitude be above a minimum, taken conservatively to be
$d_{\mathrm{8\,EeV}}>5\%$ for $E>8\,\mathrm{EeV}$. This value is around
$2.5\sigma$ below the Auger value (see Fig. 6) so it constitutes approximately
a $99\%$ C.L. lower limit on the dipole amplitude. This requirement is easily
fulfilled by many samples, especially for lower densities.

Table 1 reports the number of cases, out of the 1000 realizations for each
source density, fulfilling the two criteria. The combination of both criteria
excludes densities $n\leq 10^{-3.5}\,{\rm Mpc}^{-3}$ at more than 99% C.L.,
for the case of negligible EGMF being treated in this section.

$n/{\rm Mpc}^{-3}$ | $10^{-5}$ | $10^{-4}$ | $10^{-3.5}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$
---|---|---|---|---|---|---
$d_{\mathrm{8\,EeV}}$ >5$\%$ | 906 | 802 | 717 | 587 | 416 | 288
$C_{l>1}$ in $99\%$ iso | 0 | 0 | 3 | 12 | 273 | 569
both | 0 | 0 | 0 | 5 | 100 | 152
Table 1: The entries in the table give the number of realizations out of 1000,
which satisfy the criterion in the first column, for the
$\beta_{\mathrm{EGMF}}=0$ case. "$C_{l>1}\,99\%$ iso" means that apart from
the dipole, all multipoles are within the 99% isotropic expectation for all
energy bins and thresholds, as observed for the data (c.f., Fig. 10).

Figure 10: Angular power spectrum in the baseline model, for different source
number densities $n$, for $E>8\,\text{EeV}$ (top) and $E>32\,\text{EeV}$
(bottom). The markers show the median value for the given density and the
error bars indicate the range encompassing $68\%$ of the variations from
statistical fluctuations in source distribution and cosmic ray samples. The
99% expectation from isotropic arrival directions is shown as a red line, with
the mean $\pm$ standard deviation shown as a grey band, as in Caccianiga, L.
for the Pierre Auger and Telescope Array Collaborations (2023). Note that for
Auger and combined Auger+TA data, all $C_{l>1}$ are within isotropic
expectations, i.e., below the red line, apart from the dipole. The dipole
amplitudes for different densities, and for the data, are shown in Fig. 6 in
the Appendix.

For reference, the power spectrum predicted in the baseline model is displayed
in Fig. 10 for the two cumulative energy bins, $>$8 EeV and $>$32 EeV, and
four source number densities. For good comparability, the same style is used
to mark the expectations from isotropic simulations as used by the combined
working group between Auger and TA (Caccianiga, L. for the Pierre Auger and
Telescope Array Collaborations, 2023). One sees from Fig. 10 that the level of
anisotropy in the arrival flux increases significantly with decreasing source
number density. This is as expected: the smaller the number of sources, the
greater the impact of the few local sources on the total flux.

The power spectrum for the JF12-reg model is displayed in Appendix C. There
are sizeable quadrupole and octupole moments for all densities. That means
that if the magnetic field is solely coherent without any turbulence from GMF
or EGMF, it is not possible to obtain a sizeable dipole moment without
significant higher multipole moments. In order to reproduce the data, some
level of turbulence, from either GMF or EGMF, is hence absolutely necessary.

###  5.3 Combined constraints on the EGMF
and source number density

As was demonstrated in the previous subsection, the source number density has
a strong impact on the anisotropy level of UHECRs, specifically the power
spectrum and dipole amplitude. These observables are also impacted by random
fields along the UHECRs’ trajectories as we saw in Sec. 5.1. Therefore, in
this section, we perform a combined analysis considering a range of values for
$\beta_{\mathrm{EGMF}}=B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$ as well as
the source number density. We again create 1000 realizations of each scenario
based on the baseline and JF12-reg models and determine how many of them
fulfill the two criteria outlined in Sec. 5.2: sufficient dipole amplitude but
higher multipoles not exceeding the 99% isotropic C.L.. It is important to
note that the dipole amplitude is affected less by turbulent magnetic fields
than higher multipoles (de Oliveira et al., 2023), so that we do not expect a
linear dependency between the allowed values of $\beta_{\mathrm{EGMF}}$ and
$n$.

Figure 11: Combined constraints on the the source number density $n$ and EGMF
parameter $\beta_{\mathrm{EGMF}}\equiv
B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$; dark grey regions are favored. The
intensity bar at the top shows the number of simulations out of 1000 total for
the baseline model that have both a sufficiently large dipole and higher
multipole moments small enough to be compatible with the 99% isotropic
expectations as found for the data. The white region is excluded at more than
99% C.L.. The dashed region indicates the excluded region with no Galactic
random field. Two extreme cases discussed in the text are depicted by blue
markers. Characteristic estimates of the number densities of some steady
source candidates are shown with red lines, and indicative locii of transient
source candidates with orange lines; see Sec. 5.4 for discussion and
references.

In Fig. 11, the number of simulations that fulfill both criteria combined is
depicted graphically. The white area is excluded at more than 99% C.L. as
there are no simulations out of 1000 that fulfill both criteria at the same
time. The lower left part of the phase space is excluded by too-large
multipole moments, and the upper right part by a too-small dipole amplitude. A
clear correlation of the allowed area is visible: if the EGMF is small with
$\beta_{\mathrm{EGMF}}\lesssim 0.1$, the source density has to be large with
$n\gtrsim 10^{-4}\,{\rm Mpc}^{-3}$. This part of the phase space also includes
the case of no EGMF discussed in Sec. 5.2, indicated with a circle marker in
Fig. 11. If, instead, the EGMF is sizeable – $\beta_{\mathrm{EGMF}}\gtrsim
0.3$ – both large source densities $n\gtrsim 10^{-3}\,{\rm Mpc}^{-3}$ and very
small densities $n\lesssim 10^{-9}\,{\rm Mpc}^{-3}$ are excluded at 99% C.L.,
producing too small dipole amplitude or too large multipoles in the respective
cases. An EGMF stronger than $\beta_{\mathrm{EGMF}}\geq 3$ leads to too small
dipole moments for any source number density and is thus excluded at 99% C.L..
Note that this is a competitive upper limit on the EGMF (Durrer & Neronov,
2013). It is in the same range as the 95% C.L. upper limit of 1 nG for 1 Mpc
coherence length derived in Pshirkov et al. (2016). As the EGMF
$\beta_{\mathrm{EGMF}}$ and/or UHECR source density become better constrained
by other means, the constraints from UHECR anisotropies embodied in Fig. 11
will further restrict the other.

Another scenario that can be excluded based on the combination of sufficient
dipole magnitude and absence of higher multipoles is a homogeneous source
distribution in the case of negligible EGMF. For a sufficiently low source
density, the required dipole magnitude can result from the local inhomogeneity
of sources arising randomly in an otherwise homogeneous distribution (Harari
et al., 2015). However, our analysis shows that the required small source
density leads to too-large higher multipole moments, and there is no density
for which both criteria can be fulfilled simultaneously. This is shown in Sec.
D of the Appendix. For non-negligible EGMF with $\beta_{\mathrm{EGMF}}\simeq
1$, it is in principle possible to reproduce the dipole amplitude while
keeping all higher multipoles within the 99% isotropic C.L. (see Fig. 16 in
the Appendix) but much less likely for both criteria to be satisfied at the
same time than in the LSS model. Additionally, the direction is obviously only
very rarely coincident with the data for a random fluctuation of a homogeneous
source distribution.

###  5.4 Astrophysical Implications

The obtained limits on the source number density enable conclusions on the
possible classes of UHECR sources. The densities of different classes of
objects can be measured by astronomical observations. For reference, the
number of Milky Way-like galaxies is $\mathcal{O}(10^{-2})$ per cubic Mpc
(e.g., Conselice et al. (2016)). The number density of AGN is a strong
function of their luminosity333The number density of AGN evolves strongly with
the redshift; the values quoted here are at $z=0$, as appropriate to UHECRs in
our energy range which cannot propagate over cosmological distances. See Zaw
et al. (2019) for a complete, homogeneous catalog of local AGN.. Low-
luminosity AGN are quite common, with 20-40% of all galaxies in SDSS at low
redshift being broad-line AGN or satisfying the Kauffman criterion for narrow-
line AGN (Zaw et al., 2019); see (Ho, 2008; Danforth et al., 2016) for earlier
work. At high-luminosity, the density is $\mathcal{O}(10^{-5}-10^{-6})\,{\rm
Mpc}^{-3}$ (Gruppioni et al., 2013; Best & Heckman, 2012); blazars have been
estimated to have a density as sparse as $\sim\mathcal{O}(10^{-7})\,{\rm
Mpc}^{-3}$ (Ajello et al., 2013). Starforming galaxies also containing an AGN
are very common with density $\mathcal{O}(10^{-3}-10^{-2})\,{\rm Mpc}^{-3}$
(Gruppioni et al., 2013), however the highest luminosity starburst galaxies
(also known as ULIRGs, for Ultra-Luminous Infrared Galaxies) are much rarer
with density around $\mathcal{O}(10^{-5}-10^{-4})\,{\rm Mpc}^{-3}$ (Gruppioni
et al., 2013; Murase & Fukugita, 2019; Alves Batista et al., 2019).

Thus if the EGMF smearing is insignificant, i.e.,
$\beta_{\mathrm{EGMF}}\lesssim 0.1$, our limit on the source density of $n\geq
10^{-3.5}\,{\rm Mpc}^{-3}$, rules out very sparse object classes such as high-
luminosity AGNs and ULIRGs. In this case, UHECRs must come from more common
source types such as starforming / normal galaxies or low-luminosity AGNs, or
else transient sources (see below). If, however, the EGMF produces more
smearing, so $\beta_{\mathrm{EGMF}}\gtrsim 0.3$, sources have to be sparse in
order to produce a sizeable enough dipole. With EGMF strengths in the range of
$\beta_{\mathrm{EGMF}}\simeq 0.5-1$, the sources could be e.g. starburst
galaxies, high-luminosity AGNs, or blazars. This is also in agreement with the
finding of Eichmann et al. (2022) that radio galaxies could explain the
observed dipole for $\beta_{\mathrm{EGMF}}\simeq 1$. Our constraints imply
that $\beta_{\mathrm{EGMF}}\gtrsim 3$ is highly disfavored, except in the
context of a specially crafted scenario.

However, another possibility is that UHECRs are produced in transient events,
for instance giant AGN flares, tidal disruption events (TDEs) (Farrar &
Gruzinov, 2009), collapsars (long gamma-ray bursts) (LGRBs) (Waxman & Miralda-
Escudé, 1996) or compact binary mergers (Neutron Star (NS) - Black Hole (BH)
or binary NS).

The estimated $z=0$ rate of long GRBs is
$\Gamma_{\mathrm{long\,GRBs}}=1.3\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$
(Wanderman & Piran, 2010). However, given the lack of neutrino-GRB
correlations (Abbasi et al., 2022) and also the fact that the total power of
long GRBs in UHECRs would have to exceed their gamma-ray luminosities by more
than an order of magnitude (Farrar & Gruzinov, 2009), at present long GRBs are
disfavored as UHECR sources.

Meanwhile, TDEs have become an even more attractive option thanks to
discoveries of coincidences between TDEs and high energy neutrinos (Stein et
al., 2021) and theoretical developments (Farrar & Piran, 2014; Piran &
Beniamini, 2023); AGN flares as a possible source class have received little
follow-up. And with recent evidence of jet formation in NS-BH and BNS mergers,
these are also potential accelerators with NS-BH especially attractive on
account of having less baryon loading.

Due to deflections in the EGMF during propagation, the UHECRs produced in a
given transient event arrive to Earth over an extended time period we denote
as $\tau_{\mathrm{eff}}$. Therefore the effective number density of sources
contributing at any given epoch is

$n_{\rm eff}\approx\Gamma\,\tau_{\rm eff}~{},$ (5)

where $\Gamma$ is the rate of the transient events which produce the UHECRs.
From (Achterberg et al., 1999) the time-spreading is
$\mathcal{D}_{0}D^{2}/(3c)$

$\tau_{\rm
eff}=0.14\left(\frac{D}{\mathrm{Mpc}}\frac{\mathrm{EV}}{R}\beta_{\mathrm{EGMF}}\right)^{2}\\!\mathrm{Myr}\\\
=34\,\beta_{\mathrm{EGMF}}^{2}\,\mathrm{Myr}.$ (6)

In the last equality we evaluated the expression for rigidity
$R=4.5\,\mathrm{EV}$ and distance $D\simeq 70$ Mpc (see Fig. 1 and 4),
appropriate to the more constraining $E\geq 32$ EeV dataset.

The total rate of TDEs is estimated to be
$\Gamma_{\mathrm{all\,TDE}}=10^{2-3}\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$
(van Velzen & Farrar, 2014; Teboul et al., 2023; Yao et al., 2023), but not
all TDEs may be capable of producing UHECRs. The rate of jetted TDEs is
difficult to determine robustly because the radio emission signaling the
presence of a jet may only emerge many years after the event (Piran &
Beniamini, 2023). With this in mind the estimate of
$\Gamma_{\mathrm{jetted\,TDE}}=0.03\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$
(Andreoni et al., 2022) can be considered a lower bound on the jetted TDE
rate. The LIGO-Virgo-Kagra collaboration estimates the rate of NS-BH mergers
as (7.8-140)$\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$, with the BNS rate
possibly greater than an order of magnitude larger (Abbott et al., 2023).

In Fig. 11, we show with orange lines the relation of $\beta_{\mathrm{EGMF}}$
and
$n_{\mathrm{eff}}=\Gamma\tau_{\mathrm{eff}}=34\,\Gamma\beta_{\mathrm{EGMF}}^{2}$
following Eq. 6 for the different transient candidates discussed above and
using $30\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$ as a representative choice for
NS-BH mergers. It is visible that all of them cross the allowed region, so
none can be excluded as sources from this study. In general, transient source
candidates with the lowest rates require stronger EGMF time-averaging in order
to have a sufficient number density of sources contributing at any given
epoch.

##  6 Summary and Conclusions

In this work we explored the implications of the hypothesis that UHECR sources
follow the Large Scale Structure of the local universe, improving upon the
earlier study by Ding et al. (2021) by implementing a self-consistent
treatment of propagation effects and a simultaneous fit to the UHECR spectrum,
composition, and dipole data. The best-fit parameters of the source injection
are in agreement with previous works using a homogeneous source distribution
(The Pierre Auger Collaboration, 2023a, b), with a hard spectral index and
predominantly nitrogen injection. We verified that the results are almost
completely independent of systematic uncertainties from the hadronic
interaction model and variations of the source injection. We investigated the
sensitivity of various features of the data to uncertainties in the strength
of the Galactic random field.
$\bullet$ We demonstrated that the _dipole amplitude_ and its energy evolution
are well reproduced by the LSS model. The measured amplitude is bracketed by
our two tested GMF models, chosen to encompass the uncertainty in the random
component of the GMF: the baseline JF12 model (Jansson & Farrar, 2012b, a)
including a turbulent component with 30 kpc coherence length, and “JF12-reg"
with only the coherent field part.
$\bullet$ We demonstrated that the _dipole direction_ is independent of the
turbulent part of the GMF. It is, however, strongly influenced by the regular
component of the GMF and not perfectly reproduced using JF12. In the future,
the dipole direction can help select among different coherent magnetic field
models (Unger & Farrar, 2023), largely independently of the turbulent field
part.

We investigated the possibility of a _bias_ between UHECR sources and the LSS.
Based on the ability to reproduce the observed dipole, we find that the LSS
matter distribution is a good, bias-free proxy for the UHECR distribution. In
particular, UHECR sources reside in both high- and intermediate-density
regions; we came to no definite conclusion regarding their presence in low-
density regions.

We used the amplitude of the dipole and upper limits on higher multipoles to
constrain the extragalactic magnetic field and the number density of UHECR
sources. The arrival direction smearing due to the EGMF depends on the
combination of rms field strength, $B$, and its coherence length $L_{c}$,
through the combination $\beta_{\mathrm{EGMF}}\equiv B\sqrt{L_{c}}$.
$\bullet$ In the “continuum limit" of many low-luminosity UHECR sources, we
place the 99% C.L. upper limit $\beta_{\mathrm{EGMF}}\lesssim
0.2\,\,\mathrm{nG\,Mpc^{1/2}}$. If the extragalactic smearing were larger than
this, the strength of the dipole amplitude could not be as large as observed,
within the LSS framework.
$\bullet$ If the arrival direction smearing by the EGMF is negligible, the
number density of UHECR sources cannot be too small: $n>10^{-3.5}\,{\rm
Mpc}^{-3}$ at 99% C.L.. If the sources are more rare than this, the power
spectrum would have significant quadrupole and octupole moments, contrary to
observation. Compared to the previous best limit of
$n\gtrsim\mathcal{O}(10^{-4})\,{\rm Mpc}^{-3}$ (The Pierre Auger
Collaboration, 2013), our limit is not only more constraining but also more
accurate in that it takes into account deflections in the Galactic magnetic
field and consistently treats the mixed composition. The greatly enhanced
event statistics accumulated in the decade since that work call for use of the
more realistic LSS source distribution as done here.
$\bullet$ We set combined constraints in the parameter space of the source
density $n$ and the EGMF smearing parameter $\beta_{\mathrm{EGMF}}$, shown in
Fig. 11. A correlation between the two quantities is visible.
$\beta_{\mathrm{EGMF}}$ up to $\approx 1\sqrt{\mathrm{Mpc}}\ \mathrm{nG}$ is
possible, but requires the sources to be rare with density $n$ in the range
$\approx(10^{-9}-10^{-4})\,{\rm Mpc}^{-3}$.
$\bullet$ We investigated the idea that the dipole can arise as a statistical
fluctuation of a homogeneous source distribution (Harari et al., 2015). In
general that scenario is in conflict with the absence of higher multipoles,
but it can work for a sufficiently rare source population $n\lesssim
10^{-5}\,{\rm Mpc}^{-3}$ and relatively strong EGMF
$\beta_{\mathrm{EGMF}}\simeq(0.1-2)$, albeit in a fairly small fraction of
trials (see Fig. 16). The direction of the dipole is of course totally random
and not predictable or interpretable in that scenario, unlike for the LSS
model.
$\bullet$ An EGMF smearing parameter $\beta_{\mathrm{EGMF}}\geq
3\sqrt{\mathrm{Mpc}}\ \mathrm{nG}$ is excluded at 99% C.L. for all source
densities, for both LSS and homogeneous sources. This places a competitive
upper limit on the EGMF (Durrer & Neronov, 2013) which will be investigated in
more depth in the future.

We placed limits on the effective source density, which restrict the possible
classes of UHECR sources. If the EGMF smearing is small, UHECR sources must be
common – for example low-luminosity AGNs, star forming galaxies, or normal
galaxies. Rare, powerful sources like high-luminosity AGNs, blazars, or
luminous starburst galaxies can be excluded as the exclusive source of UHECRs
unless there is sufficient EGMF smearing (but that is limited by the magnitude
of the dipole – see above). Transient sources can also be constrained by our
bounds, if their rate is accurately enough known, since their effective number
density increases as the square of the spreading parameter
$\beta_{\mathrm{EGMF}}$ times the rate.

Finally, we investigated correlated anisotropies in flux, mean rigidity, and
composition (see Appendix). The effects are too small to be statistically
significant at present, but may be detectable in the not-too-distant future
thanks to ongoing data taking and upgrades for better mass determination by
Auger and TA.

## Acknowledgements

We are grateful for fruitful discussions with many colleagues including Eric
Mayotte, Marco Muzio, Foteini Oikonomou, Michael Unger, and members of the
Pierre Auger Collaboration. Also, we thank Chen Ding for sharing insights into
his work, Noemie Globus for sharing the CosmicFlows 2 density distribution,
and Esteban Roulet and Silvia Mollerach for comments on the manuscript. The
work of T.B. is supported by a Radboud Excellence fellowship from Radboud
University in Nijmegen, Netherlands and that of G.R.F. is supported by
National Science Foundation Grant No. PHY-2013199. T.B. thanks the Center for
Cosmology and Particle Physics of New York University for its kind
hospitality, including support from the Alfred P. Sloan Foundation
facililtating this research and the production of this paper.

##  Appendix A Influence of model uncertainties on the results

In this section, the influence of the most important model uncertainties on
the fit results is investigated.

###  A.1 Fitted model parameters

The fitted parameters for the different tested models are given in Table 2.

| baseline | JF12 reg | homogeneous |  no $\ln\mathcal{L}_{d}$ | Sibyll | sec1 | sec2 | sec3
---|---|---|---|---|---|---|---|---
sources | LSS | LSS | homogenous | LSS | LSS | LSS | LSS | LSS
HIM | EPOS-LHC | EPOS-LHC | EPOS-LHC | EPOS-LHC | Sibyll2.3d | EPOS-LHC | EPOS-LHC | EPOS-LHC
cutoff $f_{\mathrm{cut}}$ | b.e. | b.e. | b.e. | b.e. | b.e. |  sec($x$) |  sec($x^{2}$) |  sec($x^{3}$)
GMF | reg+rand | reg only | reg+rand | reg+rand | reg+rand | reg+rand | reg+rand | reg+rand
$\gamma$ | $-1.17$ | $-1.23$ | $-1.34$ | $-1.12$ | $-0.73$ | $-1.10$ | $0.74$ | $1.08$
$\log_{10}R_{\mathrm{cut}}$ | 18.2 | 18.2 | 18.2 | 18.2 | 18.2 | 18.2 | 18.8 | 18.9
$I_{\mathrm{H}}$ | 0.01 | 0.02 | 0.02 | 0.01 | 0.00 | 0.00 | 0.06 | 0.07
$I_{\mathrm{He}}$ | 0.27 | 0.27 | 0.24 | 0.27 | 0.27 | 0.28 | 0.11 | 0.00
$I_{\mathrm{N}}$ | 0.57 | 0.56 | 0.58 | 0.57 | 0.55 | 0.57 | 0.69 | 0.80
$I_{\mathrm{Si}}$ | 0.12 | 0.12 | 0.13 | 0.12 | 0.14 | 0.11 | 0.12 | 0.07
$I_{\mathrm{Fe}}$ | 0.01 | 0.04 | 0.04 | 0.03 | 0.04 | 0.03 | 0.03 | 0.06
$\nu_{X_{\mathrm{max}}}/\sigma$ | $-0.88$ | $-0.93$ | $-0.95$ | $-0.91$ | $-0.84$ | $-0.84$ | $-1.26$ | $-1.50$
$\ln\mathcal{L}_{E}$ | $-92.3$ | $-91.1$ | $-92.4$ | $-92.3$ | $-91.4$ | $-92.2$ | $-94.2$ | $-99.5$
$\ln\mathcal{L}_{X_{\mathrm{max}}}$ | $-228.7$ | $-229.2$ | $-229.3$ | $-228.7$ | $-229.1$ | $-229.0$ | $-230.4$ | $-230.0$
$\ln\mathcal{L}_{d}$ | 12.1 | $11.8$ | $-8.5$ |  12.0* | 11.2 | 12.1 | 10.6 | $9.5$
$\ln\mathcal{L}_{\mathrm{syst}}$ | $-0.4$ | $-0.4$ | $-0.5$ | $-0.4$ | $-0.4$ | $-0.4$ | $-0.8$ | $-1.1$
$\ln\mathcal{L}_{\mathrm{sum}}$ | $-309.3$ | $-308.9$ | $-330.6$ | $-321.4$ | $-309.6$ | $-309.4$ | $-314.9$ | $-321.1$
Table 2: Best-fit source parameters and likelihood values for the multiple
tested models. The differences between models are indicated in the top rows
(see also text), with HIM relating to the hadronic interaction model and b.e.
to the broken exponential cutoff. The baseline model results and
interpretation are presented in Sec. 3. Some of the other tested models are
discussed in Appendix A. (* not included in $\ln\mathcal{L}_{\mathrm{sum}}$)

###  A.2 Uncertainty from the air shower model

To test the influence of the hadronic interaction model on the result, we
compare a model using Sibyll2.3d (Riehn et al., 2020) to the baseline one
using EPOS-LHC. The best-fit parameter values for that case are given in Table
2. The fit parameters as well as the modeled observables (see Fig. 12), are
almost independent of the hadronic interaction model. Only the spectral index
at injection is softer in the case of Sibyll, as was also found in The Pierre
Auger Collaboration (2023a). The likelihood value of the Sibyll model is
almost the same as the baseline one. This is not as found in The Pierre Auger
Collaboration (2023a), where the fit using Sibyll is significantly worse.
Note, however, that compared to (The Pierre Auger Collaboration, 2023a) we use
CRPropa instead of Simprop, and also a higher energy threshold of 8 EeV.

Figure 12: Best-fit energy spectrum (left), first two moments of the
$X_{\mathrm{max}}$ distributions (middle), and dipole amplitudes (right) for
the Sibyll model (see best-fit parameters in Table 2), using Sibyll 2.3d as
hadronic interaction model.

As the dipole amplitude and direction are almost not influenced by the
hadronic interaction model, the conclusions in Sec. 5, Sec. 5.2, and Appendix
F can be made independently of these uncertainties.

###  A.3 Influence of the cutoff function shape at injection

To investigate the influence of the choice of cutoff function on the fit
result, we vary the functional form of $f_{\mathrm{cut}}$ as described in Sec.
2.1. The resulting fit parameters are given in Table 2. This is important
because of the small baseline best-fit maximum source rigidity of
$R_{\mathrm{cut}}=10^{18.2}$ V - hence the cutoff of the dominant nitrogen
contribution already happens at $E_{\mathrm{cut}}^{N}=7\times
10^{18.2}\,\text{eV}\approx 10\,\text{EeV}$. This means that the source
injection in the energy range $E>8$ EeV studied in this work is dominated by
the cutoff region instead of the power-law with spectral index $\gamma$.

For the soft cutoff with $f_{\mathrm{cut}}(x)=\text{sec}(x)$, almost no
difference can be seen in the fit parameters and likelihood values compared to
the broken exponential cutoff. When the steeper cutoff functions with
$\text{sec}(x^{2})$ and $\text{sec}(x^{3})$ are used, the injected spectrum
becomes significantly softer, leading to positive values of $\gamma$ more in
agreement with expectations of shock acceleration. The same is also observed
in Gonzales, J.M. Gonzales for the Pierre Auger Collaboration (2023). The
softening of the injected spectrum happens because the pronounced features of
the measured spectrum at Earth are now modeled by the steep cutoff instead of
a hard injection spectrum. As the spectral index is correlated with the
maximum rigidity, this also leads to higher values of $R_{\mathrm{cut}}$. The
integral element contributions are relatively independent of the cutoff shape.

The value of the log-likelihood decreases for the steeper cutoffs, indicating
a worse description of the data. This is mainly due to the energy spectrum,
which cannot be described as well with the softer injected spectra.
Additionally, the dipole likelihood $\mathcal{L}_{d}$ decreases with the
steepening of the injection cutoff. The reason for this is the larger mixing
of elements induced by the softer spectral index. This leads to stronger
contributions of heavier elements at smaller energies as in the baseline case,
which diffuse stronger. This induces a further decrease of the amplitude of
the modeled dipole in the lower energy bins compared to the baseline model.
Based on these results, we only use the broken exponential injection for the
main part of this work.

##  Appendix B Source number density influence on dipole amplitude and
direction

Here we investigate how well the observed dipole amplitude and direction are
reproduced for different source number densities for the baseline GMF model,
without EGMF, using the same setup as in Sec. 5.2. The variation of the dipole
direction and amplitude are shown in Fig. 13. For small densities, $n\lesssim
10^{-4}\,{\rm Mpc}^{-3}$, only a few of the sparse sources make significant
contributions to the total flux at Earth, so that the variations between
realizations are sizeable. This leads to a big spread of dipole directions and
generally too-large dipole amplitudes. For large densities, $n\gtrsim
10^{-3}\text{Mpc}^{-3}$, the variations due to the number of events dominate
over the variations in the source distribution, so that the dipole amplitude
and directions mostly agree with the continuous case.

Figure 13: Spread of dipole direction (upper row) and dipole amplitude $A$
(lower row) for various source number densities $n$, for two different energy
thresholds. The pink contours in the skymaps and filled region in the
histogram indicate the statistical variations due to the limited number of
events for the continuous (baseline) model, as in Fig. 5.

Because even in the continuous case, the agreement between the predicted and
measured dipole direction is not perfect (especially at lower energies, see
Fig. 5), it is currently not useful to characterize the source density by
measuring the number of realizations that reproduce the observed dipole
direction. So, we choose to outline the intended method that could be used in
future works with updated GMF, where the continuous model gives an accurate
description of the dipole at all energies: instead of characterizing the
fraction inside the measured 68% region, here, the fraction inside the modeled
statistical uncertainty will be used.

Figure 14: Fraction of realizations where dipole amplitude (left) / direction
(middle) / both at the same time (right) are within the 68% uncertainty of the
continuous model depending on the source number density $n$. Three separate
energy bins are shown, as well as a red line indicating the combination of all
three at the same time. The black line at 68% indicates the expected value for
infinite density for each separate energy bin.

In Fig. 14, the fraction of realizations where the dipole amplitude and/or
direction is within the $68\%$ statistical uncertainty of the continuous model
(from the measured number of events, same as Fig. 5 and 6) is depicted. It is
evident that for small densities $n\lesssim 10^{-4}\,{\rm Mpc}^{-3}$, the
amplitudes become too large, and the dipole directions too random as was also
visible in Fig. 13. For $n=10^{-4}\,{\rm Mpc}^{-3}$, only 1 out of the 1000
realizations reproduces both measures at the same time for all energy bins.
Even smaller densities are excluded completely, mostly due to the too-strong
dipole amplitudes.

##  Appendix C Power spectra without turbulent GMF component

In Fig. 15, the power spectrum for the JF12-reg model without turbulent GMF
component is depicted for different sampled number densities $n$. It is
visible how, compared to the case with turbulent GMF part (Fig. 10), the
quadrupole and octupole moments are much more pronounced. The pronounced
features are also visible in the predicted arrival directions for the
continuous JF12-reg model (see Fig. 5, fourth row). For no value of $n$ is a
simulation found that stays within the 99% C.L. isotropic expectation for any
density. Thus, we can exclude the case of no turbulence (neither from GMF nor
EGMF) at the 99% C.L..

Figure 15: Power spectrum for different number densities $n$; same as Fig. 10,
but for the JF12-reg model without turbulent GMF component.

##  Appendix D Alternative to LSS? A homogeneous source distribution

In The Pierre Auger Collaboration (2018), based on the results of Harari et
al. (2015), it was suggested that homogeneously distributed sources (not
following the LSS) with a density around $n\simeq 10^{-4}\,{\rm Mpc}^{-3}$,
could also explain the measured dipole energy evolution. We investigate this
by changing the source distribution to a fully homogeneous setup and repeating
all steps described in Sec. 5. The best-fit composition and spectrum at the
sources in the homogeneous case are similar to the baseline ones, see Table 2.

Analogously to Fig. 11, we show in Fig. 16 the number of samples that
simultaneously fulfill the two criteria explained in Sec. 5.2 (dipole
amplitude above 5% for $E>8$ EeV and all higher multipole moments inside 99%
isotropic C.L.) for a homogeneous underlying source distribution. As the
excluded regions for the models with and without Galactic random field are
similar, all following conclusions apply independently of the Galactic
turbulence.

Figure 16: Combined constraints on $\beta_{\mathrm{EGMF}}$ and the number
density $n$; same as Fig. 11, but for a homogeneous source distribution not
following the LSS. Again, the dashed region stands for the non-excluded part
in the case without Galactic random field while the color scale shows the
number of samples for the baseline model with Galactic random field. This
shows that it is possible to explain the magnitude of the dipole and absence
of higher multipoles as a statistical fluctuation of homogeneous sources, for
a sufficiently rare source population and relatively strong EGMF, albeit in a
fairly small fraction of trials (see Fig. 11 for comparison to the LSS model).
However, the direction of the dipole is totally random and not predictable or
interpretable in this scenario, unlike for the LSS model.

In the case of negligible EGMF – $\beta_{\mathrm{EGMF}}\lesssim 0.05$ – no
samples are found that fulfill both criteria. This is because relatively small
densities of $n\lesssim 10^{-3}\,{\rm Mpc}^{-3}$ are necessary to produce a
sizeable dipole from a homogeneous source distribution. But, for these small
densities, the higher multipole moments are not compatible with isotropy. This
can also be seen in the power spectra in Fig. 17.

Figure 17: Power spectrum for different number densities $n$; same as Fig. 10,
but for a homogeneous source distribution not following the LSS. The EGMF is
zero, and the GMF is JF12 reg+rand.

If the EGMF is not negligible – $\beta_{\mathrm{EGMF}}\gtrsim 0.1$ – it is
possible to fulfill both criteria at the same time (Note that, as in the LSS
case, strong EGMFs with $\beta_{\mathrm{EGMF}}\gtrsim 3$ are excluded). It is
however significantly less likely than in the LSS case (compare Fig. 16 to
Fig. 11) as it is harder to produce a sizeable dipole from a homogeneous
source distribution than from the LSS. The highest number of samples that
fulfill both criteria simultaneously is 64 out of 1000 for the homogeneous
source model, compared to $200/1000$ for the LSS model.

In addition, the dipole direction is roughly uniformly distributed in the
homogeneous source case, which makes it much less likely to align with the
measured one. This is visible in Fig. 18. Note that the dipole directions are
not completely uniformly distributed due to the anisotropic extragalactic
source distribution from sampling, in combination with preferred deflections
of the GMF.

Figure 18: Dipole directions of 100 samples for different number densities
$n$; same as Fig. 13 (upper), but for a homogeneous source distribution not
following the LSS. The EGMF is zero, and the GMF is JF12 reg+rand.

As for the homogeneous model, compared to the LSS one, significantly smaller
source densities are necessary, it becomes unlikely that dense source
candidates like Milky-Way-like / starforming galaxies or low-luminosity AGNs
are the sources of UHECRs in that case. Also, the frequency of all TDEs is too
large to be compatible. Even for the extreme jetted case, the likelihood of
being compatible with data for a homogeneous underlying distribution is small,
same as for long GRBs. Thus, in the case of homogeneous source distribution,
sparse or infrequent source candidates like blazars are favored. This is
different to the LSS model discussed in sec. 5.4, for which dense or frequent
objects are best compatible with the dipole and higher multipoles.

##  Appendix E EGMF and source number density constraints

In table 3, the number of simulations (out of 1000) that fulfill the two
criteria stated in the second column of the table, are given for different
values of $\beta_{\mathrm{EGMF}}$. For more details see sec. 5.2 and sec. 5.3.

| density $n/{\rm Mpc}^{-3}$ | $10^{-8}$ | $10^{-7}$ | $10^{-6}$ | $10^{-5}$ | $10^{-4}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$
---|---|---|---|---|---|---|---|---|---
$\beta_{\mathrm{EGMF}}=0$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 995 | 998 | 989 | 906 | 802 | 587 | 416 | 288
|  $C_{l>1}$ in $99\%$ isotropic expectation | 0 | 0 | 0 | 0 | 0 | 12 | 273 | 569
| both combined | 0 | 0 | 0 | 0 | 0 | 5 | 100 | 152
$\beta_{\mathrm{EGMF}}=0.03$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 993 | 998 | 981 | 882 | 789 | 560 | 411 | 293
|  $C_{l>1}$ in $99\%$ isotropic expectation | 0 | 0 | 0 | 0 | 0 | 36 | 344 | 635
| both combined | 0 | 0 | 0 | 0 | 0 | 20 | 137 | 176
$\beta_{\mathrm{EGMF}}=0.1$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 989 | 993 | 971 | 836 | 701 | 455 | 293 | 185
|  $C_{l>1}$ in $99\%$ isotropic expectation | 0 | 0 | 0 | 1 | 15 | 180 | 547 | 744
| both combined | 0 | 0 | 0 | 1 | 5 | 66 | 132 | 114
$\beta_{\mathrm{EGMF}}=0.3$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 958 | 959 | 859 | 586 | 394 | 172 | 50 | 9
|  $C_{l>1}$ in $99\%$ isotropic expectation | 0 | 0 | 6 | 107 | 287 | 634 | 833 | 944
| both combined | 0 | 0 | 4 | 18 | 39 | 60 | 17 | 8
$\beta_{\mathrm{EGMF}}=1$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 156 | 132 | 100 | 42 | 5 | 5 | 0 | 0
|  $C_{l>1}$ in $99\%$ isotropic expectation | 101 | 250 | 606 | 867 | 959 | 985 | 998 | 1000
| both combined | 59 | 50 | 36 | 15 | 1 | 0 | 0 | 0
$\beta_{\mathrm{EGMF}}=3$ |  $d_{\mathrm{8\,EeV}}$ >5$\%$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
|  $C_{l>1}$ in $99\%$ isotropic expectation | 186 | 343 | 720 | 959 | 985 | 994 | 995 | 998
| both combined | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Table 3: Same as table 1, for different values of $\beta_{\mathrm{EGMF}}\equiv
B/\mathrm{nG}\,\sqrt{L_{c}/\mathrm{Mpc}}$.

##  Appendix F Additional model predictions

In the following, we will discuss predictions following from the LSS model,
and what these would mean for UHECR measurements. The arrival-direction-
dependent composition anisotropy presented in Sec. F.1 has already been
investigated on data. The energy-flux correlation and the rigidity dependency
of the dipole discussed in sec. F.2 respective sec F.3 could be used for
future checks of the model’s correctness.

###  F.1 Directional composition anisotropy

In the Auger data (Mayotte, E. for the Pierre Auger Collaboration, 2021, 2023)
there is an indication for a composition anisotropy, where the lighter events
arrive predominantly in the Galactic North and South, while the heavier events
arrive closer to the Galactic plane. To investigate, if this anisotropy could
be an effect of the Galactic magnetic field, e.g. from the larger field
strengths closer to the Galactic plane, we investigate the size of possible
composition anisotropies over the sphere for the LSS model.

The normalized shower depth (Mayotte, E. for the Pierre Auger Collaboration,
2021) (corrected for the energy evolution effect) as a function of the
direction, is displayed in Fig. 19 (left). For the baseline continuous model,
the composition does not vary greatly over the sky, with a maximum difference
of $\Delta X\lesssim 2$ g/cm2. The lighter part (red) is correlated with the
flux overdensity seen in Fig. 5, which is expected since the lighter (low $Z$)
component diffuses less in the GMF. The correlation of the heavier part of the
events with the Galactic plane region, as observed in Auger data, is not seen
for the LSS model. From this, it can be concluded that it is not a generic
feature induced by the GMF, but instead, if not just a coincidence, must be
due to specific source directions. This could be investigated in a future
publication.

Figure 19: Left: Normalized (corrected for the energy evolution) maximum
shower depth $X_{\mathrm{max}}^{\mathrm{norm}}$ for the continuous baseline
model, with a $30^{\circ}$ tophat smoothing; in analogy to Mayotte, E. for the
Pierre Auger Collaboration (2021). Right: Histogram of the maximum difference
$\Delta X$ between smallest and largest $X_{\mathrm{max}}^{\mathrm{norm}}$ in
$30^{\circ}$ tophat smoothed maps such as the one on the left, as a function
of the source number density $n$, assuming 4000 events.

The effect of the source number density on the composition anisotropy is
visualized in Fig. 19 (right) for the case of no EGMF. The maximum difference
between the smallest and largest $X_{\mathrm{max}}^{\mathrm{norm}}$ in
$30^{\circ}$ tophat smoothed maps (like Fig. 19) of the 1000 randomly drawn
realizations is depicted. For the maps, we use a smaller event statistic of
$N_{\mathrm{events}}=4000$, as the analysis is based on data by the
fluorescence detector with a limited duty cycle. It is visible how only very
small densities $n\lesssim 10^{-7}\,{\rm Mpc}^{-3}$ can lead to differences
between on- and off-plane regions approaching the value of $9.8$ g/cm2
reported in Mayotte, E. for the Pierre Auger Collaboration (2021). To reach
the updated value of $5.1$ g/cm2 (Mayotte, E. for the Pierre Auger
Collaboration, 2023), still a density as small as $n\approx 10^{-6}\,{\rm
Mpc}^{-3}$, is necessary. As described above, such small densities lead to
very large flux anisotropies which are not in agreement with the data for
negligible EGMF.

We have also investigated the possibility of arrival-direction-dependent
composition anisotropies for the homogeneous source model (Appendix D), and
have found that the directions of possible anisotropies in the sky are more
random in that case. The size of the maximum $\Delta X$, however, does not
depend on the source model.

###  F.2 Directional energy anisotropy

We have also tested the possibility of a directional energy anisotropy using
our model. In Fig. 20, we show the mean energy per pixel (upper row) for the
illumination and the arrival maps for two energy thresholds for the baseline
model. For comparison, the respective flux maps are also shown (lower row). A
clear correlation is visible: the regions with larger flux also exhibit higher
energies which is an effect of the extragalactic propagation. Note that the
mean energy is dominated more by nearby structures than the flux, as can be
seen in comparison to the illumination maps for different distances shown in
Fig. 2.

Figure 20: The upper row shows the mean energy per pixel, for the illumination
(first column for $E>8$ EeV, third column for $E>32$ EeV) and the arrival
(second column for $E>8$ EeV, fourth column for $E>32$ EeV). For comparison,
the flux maps for the respective energy thresholds are depicted in the lower
row. A correlation between larger flux and higher energy is visible.

We have estimated roughly if the correlation between energy and flux could be
tested by Auger by using the simulated event samples described in sec. 5.2. We
have found that the expected correlation coefficient is relatively small
$\mathcal{O}(\lesssim 5\%)$ which is mostly due to limited statistics. Thus,
such a correlation could be a target for future observatories like GCOS (Alves
Batista et al., 2023).

###  F.3 Rigidity dependency of the dipole

Due to the improving sensitivity of UHECR observatories to the masses and
charges of individual events thanks to deep learning methods and detector
upgrades (Castellina, A. for the Pierre Auger Collaboration, 2019; Glombitza,
J. for the Pierre Auger Collaboration, 2023), it can be expected that effects
like a possible evolution of the dipole amplitude and direction with the
rigidity can be tested in the near future. The predicted dipole evolution with
the rigidity based on our LSS model is shown in Fig. 21. Here, it is clearly
visible how the highest rigidity events lead to a much stronger dipole
amplitude as they are less affected by diffusion in the GMF. The effect is
significantly stronger as a function of rigidity than it is seen for the
different energy bins in Fig. 5 and 6, as expected because magnetic field
deflections scale with rigidity. The events with the lowest rigidity even show
almost no dipolar anisotropy at all, while the highest rigidity ones arrive
significantly more in the Galactic North.

Testing for a rigidity evolution of the dipole, instead of just the energy
evolution would allow us to disentangle the effects of propagation and GMF
deflection, and would thus work as a proof of UHECR deflections in the GMF.
Additionally, identifying regions where the highest rigidity events are
expected could be advantageous for backtracking to identify individual source
candidates.

Figure 21: Predicted rigidity dependency of the flux $E>8\,\text{EeV}$. Each
map contains 25% of the expected events ordered in rigidity from lowest to
highest. The demarcation rigidities are $R=1.7\,\mathrm{EV}$,
$R=4.3\,\mathrm{EV}$, and $R=5.6\,\mathrm{EV}$. The star marker indicates the
direction of the dipole component.

## References

  * Abbasi et al. (2022) Abbasi, R., et al. 2022, Astrophys. J., 939, 116, doi: 10.3847/1538-4357/ac9785
  * Abbott et al. (2023) Abbott, R., et al. 2023, Phys. Rev. X, 13, 041039, doi: 10.1103/PhysRevX.13.041039
  * Achterberg et al. (1999) Achterberg, A., Gallant, Y. A., Norman, C. A., & Melrose, D. B. 1999.  https://arxiv.org/abs/astro-ph/9907060
  * Aharonian et al. (2023) Aharonian, F., Aschersleben, J., Backes, M., et al. 2023, ApJL, 950, L16, doi: 10.3847/2041-8213/acd777
  * Ajello et al. (2013) Ajello, M., Romani, R. W., Gasparrini, D., et al. 2013, ApJ, 780, 73, doi: 10.1088/0004-637x/780/1/73
  * Allard et al. (2022) Allard, D., Aublin, J., Baret, B., & Parizot, E. 2022, A&A, 664, A120, doi: 10.1051/0004-6361/202142491
  * Alves Batista et al. (2016) Alves Batista, R., Dundovic, A., Erdmann, M., et al. 2016, JCAP, 2016, 038, doi: 10.1088/1475-7516/2016/05/038
  * Alves Batista et al. (2019) Alves Batista, R., Biteau, J., Bustamante, M., et al. 2019, Frontiers in Astronomy and Space Sciences, 6, 23, doi: 10.3389/fspas.2019.00023
  * Alves Batista et al. (2023) Alves Batista, R., Ahlers, M., Assis, P., et al. 2023, in PoS(ICRC2023), Vol. 444, 281, doi: 10.22323/1.444.0281
  * Andreoni et al. (2022) Andreoni, I., Coughlin, M. W., Perley, D. A., et al. 2022, Nature, 612, 430, doi: 10.1038/s41586-022-05465-8
  * Best & Heckman (2012) Best, P. N., & Heckman, T. M. 2012, MNRAS, 421, 1569, doi: 10.1111/j.1365-2966.2012.20414.x
  * Blandford et al. (2018) Blandford, R., Simeon, P., & Globus, N. 2018, in 42nd COSPAR Scientific Assembly, Vol. 42, E1.5–61–18
  * Caccianiga, L. for the Pierre Auger and Telescope Array Collaborations (2023) Caccianiga, L. for the Pierre Auger and Telescope Array Collaborations. 2023, in PoS(ICRC2023), Vol. 444, 521, doi: 10.22323/1.444.0521
  * Castellina, A. for the Pierre Auger Collaboration (2019) Castellina, A. for the Pierre Auger Collaboration. 2019, in EPJ Web of Conferences, Vol. 210, 06002, doi: 10.1051/epjconf/201921006002
  * Condorelli et al. (2023) Condorelli, A., Biteau, J., & Adam, R. 2023, ApJ, 957, 80, doi: 10.3847/1538-4357/acfeef
  * Conselice et al. (2016) Conselice, C. J., Wilkinson, A., Duncan, K., & Mortlock, A. 2016, ApJ, 830, 83, doi: 10.3847/0004-637x/830/2/83
  * Danforth et al. (2016) Danforth, C. W., Stocke, J. T., France, K., Begelman, M. C., & Perlman, E. 2016, ApJ, 832, 76, doi: 10.3847/0004-637x/832/1/76
  * de Almeida, R. for the Pierre Auger Collaboration (2021) de Almeida, R. for the Pierre Auger Collaboration. 2021, in PoS(ICRC2021), Vol. 395, 335, doi: 10.22323/1.395.0335
  * de Oliveira et al. (2023) de Oliveira, C., Maia, L. P., & de Souza, V. 2023, Revisiting the implications of Liouville’s theorem to the anisotropy of cosmic rays.  https://arxiv.org/abs/2307.13095
  * Deligny & Salamida (2013) Deligny, O., & Salamida, F. 2013, Astroparticle Physics, 46, 40, doi: 10.1016/j.astropartphys.2013.05.003
  * Ding et al. (2019) Ding, C., Globus, N., & Farrar, G. R. 2019, in PoS(ICRC2019), Vol. 358, 243, doi: 10.22323/1.358.0243
  * Ding et al. (2021) Ding, C., Globus, N., & Farrar, G. R. 2021, ApJL, 913, L13, doi: 10.3847/2041-8213/abf11e
  * Durrer & Neronov (2013) Durrer, R., & Neronov, A. 2013, A&AR, 21, doi: 10.1007/s00159-013-0062-7
  * Eichmann et al. (2022) Eichmann, B., Kachelrieß, M., & Oikonomou, F. 2022, JCAP, 07, 006, doi: 10.1088/1475-7516/2022/07/006
  * Farrar (2014) Farrar, G. R. 2014, Comptes Rendus Physique, 15, 339, doi: 10.1016/j.crhy.2014.04.002
  * Farrar & Gruzinov (2009) Farrar, G. R., & Gruzinov, A. 2009, ApJ, 693, 329.  https://arxiv.org/abs/0802.1074
  * Farrar & Piran (2014) Farrar, G. R., & Piran, T. 2014, Tidal disruption jets as the source of Ultra-High Energy Cosmic Rays.  https://arxiv.org/abs/1411.0704
  * Gilmore et al. (2012) Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, MNRAS, 422, 3189, doi: 10.1111/j.1365-2966.2012.20841.x
  * Globus & Piran (2017) Globus, N., & Piran, T. 2017, ApJ, 850, L25, doi: 10.3847/2041-8213/aa991b
  * Globus et al. (2019) Globus, N., Piran, T., Hoffman, Y., Carlesi, E., & Pomarède, D. 2019, MNRAS, 484, 4167, doi: 10.1093/mnras/stz164
  * Glombitza, J. for the Pierre Auger Collaboration (2023) Glombitza, J. for the Pierre Auger Collaboration. 2023, in PoS(ICRC2023), Vol. 444, 278, doi: 10.22323/1.444.0278
  * Golup, G. for the Pierre Auger Collaboration (2023) Golup, G. for the Pierre Auger Collaboration. 2023, in PoS(ICRC2023), Vol. 444, 252, doi: 10.22323/1.444.0252
  * Gonzales, J.M. Gonzales for the Pierre Auger Collaboration (2023) Gonzales, J.M. Gonzales for the Pierre Auger Collaboration. 2023, in PoS(ICRC2023), Vol. 444, 288, doi: 10.22323/1.444.0288
  * Gruppioni et al. (2013) Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23, doi: 10.1093/mnras/stt308
  * Hackstein et al. (2016) Hackstein, S., Vazza, F., Brüggen, M., Sigl, G., & Dundovic, A. 2016, MNRAS, 462, 3660, doi: 10.1093/mnras/stw1903
  * Hale et al. (2017) Hale, C. L., Jarvis, M. J., Delvecchio, I., et al. 2017, MNRAS, 474, 4133, doi: 10.1093/mnras/stx2954
  * Harari et al. (2010) Harari, D., Mollerach, S., & Roulet, E. 2010, JCAP, 2010, 033–033, doi: 10.1088/1475-7516/2010/11/033
  * Harari et al. (2015) —. 2015, PRD, 92, doi: 10.1103/physrevd.92.063014
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Ho (2008) Ho, L. C. 2008, Annual Review of Astronomy and Astrophysics, 46, 475, doi: 10.1146/annurev.astro.45.051806.110546
  * Hoffman et al. (2018) Hoffman, Y., Carlesi, E., Pomarède, D., et al. 2018, Nature Astronomy, 2, 680, doi: 10.1038/s41550-018-0502-4
  * Jaffe (2019) Jaffe, T. R. 2019, Galaxies, 7, 52, doi: 10.3390/galaxies7020052
  * Jansson & Farrar (2012a) Jansson, R., & Farrar, G. R. 2012a, ApJ, 757, 14, doi: 10.1088/0004-637X/757/1/14
  * Jansson & Farrar (2012b) —. 2012b, ApJ, 761, L11, doi: 10.1088/2041-8205/761/1/L11
  * Kim, J. for the Telescope Array collaboration (2023) Kim, J. for the Telescope Array collaboration. 2023, in PoS(ICRC2023), Vol. 444, 244, doi: 10.22323/1.444.0244
  * Koning et al. (2005) Koning, A. J., Hilaire, S., & Duijvestijn, M. C. 2005, AIP Conference Proceedings, 769, 1154, doi: 10.1063/1.1945212
  * Mayotte, E. for the Pierre Auger Collaboration (2021) Mayotte, E. for the Pierre Auger Collaboration. 2021, in PoS(ICRC2021), Vol. 395, 321, doi: 10.22323/1.395.0321
  * Mayotte, E. for the Pierre Auger Collaboration (2023) Mayotte, E. for the Pierre Auger Collaboration. 2023, in PoS(ICRC2023), Vol. 444, 365, doi: 10.22323/1.444.0365
  * Mollerach & Roulet (2019) Mollerach, S., & Roulet, E. 2019, PRD, 99, doi: 10.1103/physrevd.99.103010
  * Murase & Fukugita (2019) Murase, K., & Fukugita, M. 2019, PRD, 99, doi: 10.1103/physrevd.99.063012
  * Muzio et al. (2019) Muzio, M. S., Unger, M., & Farrar, G. R. 2019, PRD, 100, doi: 10.1103/physrevd.100.103008
  * Oikonomou et al. (2015) Oikonomou, F., Kotera, K., & Abdalla, F. B. 2015, JCAP, 01, 030, doi: 10.1088/1475-7516/2015/01/030
  * Pierog et al. (2013) Pierog, T., Karpenko, I., Katzy, J. M., Yatsenko, E., & Werner, K. 2013, PRD, 92, 034906, doi: 10.1103/PhysRevC.92.034906
  * Piran & Beniamini (2023) Piran, T., & Beniamini, P. 2023, JCAP, 11, 049, doi: 10.1088/1475-7516/2023/11/049
  * Pshirkov et al. (2016) Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, PRL, 116, 191302, doi: 10.1103/PhysRevLett.116.191302
  * Riehn et al. (2020) Riehn, F., Engel, R., Fedynitch, A., Gaisser, T. K., & Stanev, T. 2020, PRD, 102, 063002, doi: 10.1103/PhysRevD.102.063002
  * Stein et al. (2021) Stein, R., et al. 2021, Nature Astron., 5, 510, doi: 10.1038/s41550-020-01295-8
  * Teboul et al. (2023) Teboul, O., Stone, N. C., & Ostriker, J. P. 2023, MNRAS, stad3301, doi: 10.1093/mnras/stad3301
  * The Pierre Auger Collaboration (2013) The Pierre Auger Collaboration. 2013, JCAP, 2013, 009, doi: 10.1088/1475-7516/2013/05/009
  * The Pierre Auger Collaboration (2014) —. 2014, PRD, 90, 122006, doi: 10.1103/PhysRevD.90.122006
  * The Pierre Auger Collaboration (2015) —. 2015, NIM A, 798, 172, doi: 10.1016/j.nima.2015.06.058
  * The Pierre Auger Collaboration (2017a) —. 2017a, Science, 357, 1266, doi: 10.1126/science.aan4338
  * The Pierre Auger Collaboration (2017b) —. 2017b, JCAP, 2017, 038, doi: 10.1088/1475-7516/2017/04/038
  * The Pierre Auger Collaboration (2018) —. 2018, ApJ, 868, 4, doi: 10.3847/1538-4357/aae689
  * The Pierre Auger Collaboration (2020) —. 2020, PRD, 102, 062005, doi: 10.1103/PhysRevD.102.062005
  * The Pierre Auger Collaboration (2023a) —. 2023a, JCAP, 2023, 024, doi: 10.1088/1475-7516/2023/05/024
  * The Pierre Auger Collaboration (2023b) —. 2023b, accepted for publication in JCAP.  https://arxiv.org/abs/2305.16693
  * The Planck collaboration (2016) The Planck collaboration. 2016, A&A, 594, A19, doi: 10.1051/0004-6361/201525821
  * The Telescope Array Collaboration (2008) The Telescope Array Collaboration. 2008, Nuclear Physics B - Proceedings Supplements, 175-176, 221, doi: https://doi.org/10.1016/j.nuclphysbps.2007.11.002
  * The Telescope Array Collaboration (2020) —. 2020, ApJL, 898, L28, doi: 10.3847/2041-8213/aba0bc
  * Tinyakov & Urban (2015) Tinyakov, P. G., & Urban, F. R. 2015, Journal of Experimental and Theoretical Physics, 120, 533, doi: 10.1134/s1063776115030231
  * Unger & Farrar (2023) Unger, M., & Farrar, G. R. 2023.  https://arxiv.org/abs/2311.12120
  * van Velzen & Farrar (2014) van Velzen, S., & Farrar, G. R. 2014, ApJ, 792, 53, doi: 10.1088/0004-637X/792/1/53
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wanderman & Piran (2010) Wanderman, D., & Piran, T. 2010, MNRAS, 406, 1944, doi: 10.1111/j.1365-2966.2010.16787.x
  * Waxman & Miralda-Escudé (1996) Waxman, E., & Miralda-Escudé , J. 1996, ApJ, 472, L89, doi: 10.1086/310367
  * Yao et al. (2023) Yao, Y., Ravi, V., Gezari, S., et al. 2023, ApJL, 955, L6, doi: 10.3847/2041-8213/acf216
  * Yushkov, A. for the Pierre Auger Collaboration (2019) Yushkov, A. for the Pierre Auger Collaboration. 2019, in PoS(ICRC2019), Vol. 358, 482, doi: 10.22323/1.358.0482
  * Zaw et al. (2019) Zaw, I., Chen, Y.-P., & Farrar, G. R. 2019, ApJ, 872, 134, doi: 10.3847/1538-4357/aaffaf
  * Zonca et al. (2019) Zonca, A., Singer, L. P., Lenz, D., et al. 2019, JOSS, 4, 1298, doi: 10.21105/joss.01298