# Simulations of ultra-high Energy Cosmic Rays in the local Universe and the origin of Cosmic Magnetic Fields

Stefan Hackstein, Franco Vazza, Marcus Brüggen, Jenny G. Sorce, Stefan Gottlöber

► To cite this version:

Stefan Hackstein, Franco Vazza, Marcus Brüggen, Jenny G. Sorce, Stefan Gottlöber. Simulations of ultra-high Energy Cosmic Rays in the local Universe and the origin of Cosmic Magnetic Fields. Monthly Notices of the Royal Astronomical Society, 2018, 475 (2), pp.2519-2529. 10.1093/mn-ras/stx3354. hal-01758092

HAL Id: hal-01758092

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

Submitted on 4 Jul 2025

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

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

<figure><img src="image_1.png" /></figure>
# Simulations of ultra-high energy cosmic rays in the local Universe and the origin of cosmic magnetic fields

S. Hackstein,<fnref n="1" />* F. Vazza,<fnref n="1" />,<fnref n="2" /> M. Brüggen,<fnref n="1" /> J. G. Sorce<fnref n="3" />,<fnref n="4" /> and S. Gottlöber<fnref n="4" />

Accepted 2017 December 27. Received 2017 November 28; in original form 2017 August 6

## ABSTRACT

We simulate the propagation of cosmic rays at ultra-high energies, $\gtrsim 10^{18}$ eV, in models of extragalactic magnetic fields in constrained simulations of the local Universe. We use constrained initial conditions with the cosmological magnetohydrodynamics code ENZO. The resulting models of the distribution of magnetic fields in the local Universe are used in the CRPROPA code to simulate the propagation of ultra-high energy cosmic rays. We investigate the impact of six different magneto-genesis scenarios, both primordial and astrophysical, on the propagation of cosmic rays over cosmological distances. Moreover, we study the influence of different source distributions around the Milky Way. Our study shows that different scenarios of magneto-genesis do not have a large impact on the anisotropy measurements of ultra-high energy cosmic rays. However, at high energies above the Greisen–Zatsepin–Kuzmin (GZK)-limit, there is anisotropy caused by the distribution of nearby sources, independent of the magnetic field model. This provides a chance to identify cosmic ray sources with future full-sky measurements and high number statistics at the highest energies. Finally, we compare our results to the dipole signal measured by the Pierre Auger Observatory. All our source models and magnetic field models could reproduce the observed dipole amplitude with a pure iron injection composition. Our results indicate that the dipole is observed due to clustering of secondary nuclei in direction of nearby sources of heavy nuclei. A light injection composition is disfavoured, since the increase in dipole angular power from 4 to 8 EeV is too slow compared to observation by the Pierre Auger Observatory.

**Key words:** MHD – relativistic processes – methods: numerical – cosmic rays – ISM: magnetic fields.

## 1 INTRODUCTION

Evidence for the existence of magnetic fields have been reported for all types of structures found throughout the Universe. Galaxies host magnetic fields with typical strengths of ~5–15 µG, which were measured using Faraday rotation and synchrotron emission up to redshift $z \sim 2–6$ (e.g. Vallée 2004; Bernet, Miniati & Lilly 2013; Beck 2016; Kim et al. 2016). The magnetic field in clusters of galaxies was found to be of the order ~µG (Feretti et al. 2012). Future radio observations will offer the chance to measure the magnetization at the outskirts of clusters and in filaments that connect them (Brown 2011; Araya-Melo et al. 2012; Vazza et al. 2015). A recent study has reported upper limits on the magnetic field strength of

~0.03 µG from the absence of a correlation between synchrotron emission and the large-scale structure (LSS, Brown et al. 2017; Vernstrom et al. 2017). Limits on the magnetic fields in voids were derived from the angular power spectrum, the bispectrum and the trispectrum of the cosmic microwave background ($B_{\text{void}} < 1$ nG, Trivedi, Subramanian & Seshadri 2014; Planck Collaboration XIX 2016), absence of evolution with redshift in Faraday rotation measures ($B_{\text{void}} < 1.7$ nG, Pshirkov, Tinyakov & Urban 2016), and the lack of secondary emission around blazar sources ($B_{\text{void}} > 10^{-7}$ nG, Neronov & Vovk 2010; Alves Batista et al. 2017)<fnref n="1" />. Magnetohydrodynamical (MHD) cosmological simulations have been used to evolve magnetic fields of primordial or other origin that are

<footnote n="1">¹Hamburger Sternwarte, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="2">²INAF, Istituto di Radioastronomia di Bologna, via Gobetti 101, I-41029 Bologna, Italy</footnote>

<footnote n="3">³Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="4">⁴Leibniz Institute for Astrophysics Potsdam, An der Sternwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote n="1">¹ See however discussion in Broderick, Chang & Pfrommer (2012) for a different view of the issue.</footnote>

<footnote>* E-mail: stefan.hackstein@hs.uni-hamburg.de</footnote>
amplified during structure formation and by additional dynamo
processes (e.g. Dolag 2006; Ryu et al. 2012). These simulations
produce models of cosmic magnetic fields (CMFs) that agree to
some extent with observations (e.g. Dolag, Bartelmann & Lesch
1999; Brüggen et al. 2005; Donnert et al. 2009). All amplification
scenarios have in common that they require a seed field, whose
structure, strength, and origin is unknown.

In this paper, we probe the possibility to learn about the origin
of CMFs using measurements of cosmic rays at ultra-high energies.
Previous studies on similar topics mainly focused on properties of
the Galactic magnetic field (Stanev 1997; Takami & Sato 2008)
or small-scale anisotropies (Harari, Mollerach & Roulet 2002a;
Yoshiguchi et al. 2003). Other works used unconstrained MHD
models to study the implications of CMFs on ultra-high energy
cosmic rays (UHECRs) astronomy (Sigl, Miniati & Ensslin 2003;
Sigl, Miniati & Enßlin 2004; Sigl, Miniati & Ensslin 2004; Das
et al. 2008; Kotera & Lemoine 2008; Hackstein et al. 2016). An
overview of UHECR studies using MHD simulations can be found
in Alves Batista et al. (2017). Analytical studies on the implica-
tion of CMFs on UHECR observations are provided in Harari,
Mollerach & Roulet (2000), Harari et al. (2002a), Harari et al.
(2002b), Tinyakov & Tkachev (2005), and Takami et al. (2012).

In previous work (Hackstein et al. 2016), we found strong vari-
ance in the observables of UHECRs induced by the position of, both,
observer and sources (also cf. e.g. Sigl et al. 2004). To reduce this
cosmic variance, it is necessary to use constrained MHD models that
resemble the local Universe, as has been done by Dolag et al. (2004).
They conclude that UHECR protons are reasonably deflected only
when they cross galaxy clusters, though they assumed a rather weak
field in voids of ≲10⁻¹¹ G. Our new work expands the early work
by Dolag et al. (2004) in a few ways: (a) we use the most recent
set of initial conditions by Sorce et al. (2016), which were derived
with more updated algorithms and observational constraints (see
Sec. 2.1); (b) we relied on a different numerical method: i.e. the grid-
MHD simulations with ENZO instead of smoothed-particle hydrody-
namics simulations, which gives us a better sampling of moderate
and low resolution regions; (c) we performed a survey of magnetic
field models, rather than assuming a single specific scenario.

Data suggest that cosmic rays are fully ionized nuclei that con-
stantly hit the Earth from outer space with energies that range over
11 orders of magnitude. At low energies (<10¹⁷ eV), the predom-
inant sources were found to be supernova remnants in our own
Galaxy, where charged particles experience Fermi acceleration in
magnetic shocks (e.g. Blasi 2013). UHECRs are less prone to the de-
flection in CMFs, thus they are not confined within their host galaxy
and presumably are of extragalactic origin. The sources of UHE-
CRs are currently unknown. If we assume the same acceleration
process as at low energies, the size of the source limits the maxi-
mum energy of emitted UHECRs. This is the famous Hillas criterion
(Hillas 1984) that limits the candidates for sources of UHECRs at
~&gt;10²⁰ eV to very few objects, namely radio galaxy lobes, clusters
of galaxies, active galactic nuclei, and gamma-ray bursts (e.g. Dova
2016). Recent works have reported signs of anisotropy in simula-
tions with pure proton composition and limited source density in
correlation with the LSS (di Matteo & Tinyakov 2017; Abreu et al.
2013), which are not observed in nature. They infer lower bounds
on the density of sources of ~10⁻⁴ Mpc⁻¹. Also, they conclude that
the UHECR flux cannot be dominated by protons. In this work, we
investigate the effect of different source distributions of UHECRs
on the observed arrival directions.

This article is organized as follows: In Section 2, we present
details on the simulation of the MHD models and of the propagation

of UHECRs. The results of these simulations are then discussed in
Section 3. Our conclusions are finally given in Section 4.

2 SIMULATION

2.1 Constrained initial conditions

Simulations that resemble the local Universe stem from particular
initial conditions. Unlike typical initial conditions that abide solely
by a cosmological prior, these initial conditions are additionally
constrained by local observational data that can be either redshift
surveys (Lavaux 2010; Heß, Kitaura & Gottlöber 2013) or radial
peculiar velocities of galaxies (Kravtsov, Klypin & Hoffman 2002;
Klypin et al. 2003; Sorce et al. 2014). We use the latter with a back-
ward (by opposition to forward, Heß et al. 2013; Jasche & Wandelt
2013; Kitaura 2013; Wang et al. 2014) technique (Bertschinger
1987; Hoffman & Ribak 1991, 1992; Ganon & Hoffman 1993;
van de Weygaert & Bertschinger 1996; Bistolas & Hoffman 1998;
Lavaux et al. 2008). The catalogue of constraints is fully described
in Tully et al. (2013) and the method to produce the constrained
initial conditions is summarized in Sorce et al. (2016). The process
involves various steps from the minimization of biases (Sorce 2015)
in the catalogue of peculiar velocities to the constrained realization
technique (Hoffman & Ribak 1991) to get the final product: the ini-
tial conditions. We work within the Planck cosmology framework
(Ω<sub>m</sub>=0.307, Ω<sub>Λ</sub>=0.693, h = 0.677, σ<sub>8</sub> = 0.829, Planck Collabo-
ration XVI 2014).

2.2 MHD-simulations

The MHD simulations performed in this paper have been produced
with the cosmological grid code ENZO that follows the dynamics of
dark matter with a particle-mesh *N*-body method and uses a variety
of shock-capturing Riemann solvers to evolve the gas component
(Bryan et al. 2014). The MHD equations were solved with the
method by Dedner et al. (2002). To keep ∇·*B* as low as possible,
it uses hyperbolic divergence cleaning. The fluxes at cell inter-
faces are reconstructed with the Piecewise Linear Method. They are
evolved using the local Lax–Friedrichs Riemann solver (Kurganov
& Tadmor 2000), with time integration using the total variation di-
minishing second-order Runge–Kutta scheme (Shu & Osher 1988).
The set of simulations was run on Piz-Daint (CSCS) and made use
of the recent implementation of the Dedner algorithm using CUDA
(Wang, Abel & Kaehler 2010).

To model the local Universe at z = 0, the MHD simulations started
at z = 60 with initial conditions described in Section 2.1. We sam-
pled a volume of (500 Mpc h⁻¹)³, with 512³ cells and dark matter
particles. We use this large volume in order to remove effects from
periodic boundary conditions in the constrained subregion of the
MHD simulation. Since the initial perturbation for baryonic matter
are not provided in the initial conditions of Section 2.1, we simply
initialize baryons to the uniform cosmological density, assuming an
initial zero-velocity field for baryons everywhere. Although more
accurate ways to couple baryons to dark matter perturbations since
the beginning are possible, this choice is irrelevant for the level of
details we are concerned here (e.g. Vazza et al. 2011). Full reso-
lution of the whole box is not necessary and costly, therefore only
the constrained innermost (200 Mpc h⁻¹)³ volume was further re-
fined by a factor 32 using adaptive-mesh refinement. The refinement
here follows the standard local overdensity criterion, doubling the
cell resolution whenever the local gas overdensity was three times
larger than the surroundings, up to a maximum of five levels of
<table><caption>**Table 1.** List of magnetic field models investigated in this paper. First column: name of the model; second column: physical module for the gas component; third column: generation of magnetic field. All models were simulated within a volume of (500 Mpc h⁻¹)<fnref n="3" />. In CRPROPA, we used the innermost (250 Mpc h⁻¹)<fnref n="3" /> with 1024³ data cells and a resolution of 245 kpc h⁻¹.</caption><thead><tr><th>Mnemonic</th><th>Gas physics</th><th>Magnetic field</th></tr></thead><tbody><tr><td>B=0</td><td>non-radiative</td><td>B<sub>0</sub> = 0</td></tr><tr><td>primordial</td><td>non-radiative</td><td>B<sub>0</sub> = 0.1 nG</td></tr><tr><td>primordial2R</td><td>non-radiative</td><td>((B<sup>2</sup>)<sup>0.5</sup>) = 1 nG, n<sub>B</sub> = -3</td></tr><tr><td>primordial3R</td><td>non-radiative</td><td>((B<sup>2</sup>)<sup>0.5</sup>) = 1 nG, n<sub>B</sub> = -4</td></tr><tr><td>astrophysical</td><td>cooling and AGN feedback</td><td>5 · 10<sup>58</sup> erg, z < 4; B<sub>0</sub> = 10<sup>-11</sup> nG</td></tr><tr><td>astrophysicalR</td><td>cooling and AGN feedback</td><td>10<sup>60</sup> erg, z < 4; B<sub>0</sub> = 10<sup>-11</sup> nG</td></tr><tr><td>astrophysicalIR</td><td>cooling and AGN feedback</td><td>10<sup>60</sup> erg to 5 · 10<sup>58</sup> erg, z < 1; B<sub>0</sub> = 10<sup>-11</sup> nG</td></tr></tbody></table>

refinement (≈31 kpc h⁻¹ per cell). The clusters that form in this volume closely resemble real local structures (e.g. the Centaurus, Virgo, Coma and Perseus clusters), within typical offsets of order ≤2–3 Mpc h⁻¹ which are however not crucial for the global studies we perform here.

The limited size of computer memory used for the simulation of UHECR propagation did not allow us to use the full volume of the CMF models obtained from MHD simulations. In order to minimize effects from periodic boundaries (see Section 2.3), we restricted the simulations in CRPROPA to the innermost (250 Mpc h⁻¹)<fnref n="3" /> volume and reduced the number of cells inside that volume to 1024³. The resulting resolution is then 245 kpc h⁻¹.

The use of constrained simulations of the local Universe is an important step forward compared to our previous work (Hackstein et al. 2016), where we found a large variance in the observed properties of UHECRs from observer to observer. Given the strong impact of ≤35 Mpc h⁻¹ sources of UHECRs, it is not guaranteed that the average over many observers is representative of what can be observed by the specific observer at Earth’s location. However, in these new runs placing our observer within the Local Group allows us to remove these uncertainties. At distances >100 Mpc h⁻¹, deflection and the increasing number of sources provide an UHECR flux almost independent of the exact position of distant sources. It is therefore sufficient to model the source distribution only within that distance.

Following a procedure similar to Hackstein et al. (2016), we ran several MHD simulations with different scenarios for the origin of CMFs. In the primordial model, we used a uniform initial magnetic field of strength 0.1 nG (comoving) along each axis at z = 60. In the primordial2R and primordial3R models, similar fields were generated by drawing the magnetic field from an analytically generated power-law distribution of magnetic fields, with two different slopes for the power spectrum, n<sub>B</sub> = −3 and −4, respectively (with P<sub>B</sub> ∝ k<sup>n<sub>B</sub></sup>), see Planck Collaboration XIX (2016) for details. We have generated a power-law spectrum distribution of the vector potential in the Fourier space for a 1024³ grid, randomly drawn from the Rayleigh distribution, and we have computed the magnetic field in real space as ⃗B = ∇ × A, ensuring ∇ · ⃗B = 0 by construction. We have assumed that the maximum coherence scale of the magnetic field is 500 Mpc h⁻¹ and that the minimum scale is the root grid resolution, and that the power-law of fluctuations follows the input P<sub>B</sub> power spectrum, similar to Bonafede et al. (2013). In both cases, the normalization of the spectrum of initial fluctuations is chosen such that ((B<sup>2</sup>)<sup>0.5</sup>) = B<sub>0</sub>, i.e. the rms magnetic field is equivalent to the uniform seeding case.

The astrophysical origin of CMFs was modelled as impulsive thermal and magnetic feedback in haloes where the physical gas number density exceeded a critical value of 10⁻<fnref n="2" /> cm⁻<fnref n="3" />. The thermal

energy is released as a couple of overpressurized outflows at random opposite directions from the halo centre. The feedback magnetic energy, assumed to be 50 per cent of the injected thermal energy, is released as dipoles around the centre.

In the astrophysical model, we assumed a release of 5 × 10⁵⁸ erg per feedback episode starting from z = 4; in the astrophysicalR model, we used instead a larger budget of 10⁶⁰ erg per event. Finally, in the astrophysicalIR model, we considered a mixed scenario, where we changed the energy budget from 10⁶⁰ to 5 × 10⁵⁸ erg per event from z = 1 to 0.

All runs with astrophysical scenarios for the emergence of extragalactic magnetic fields used equilibrium radiative gas cooling, assuming a fixed metallicity of Z = 0.3 Z<sub>☉</sub>. While the cooling is necessary to trigger the onset of cooling flows and start the cooling-feedback cycle in our haloes, the large-scale distribution of gas matter outside simulated haloes is similar across all runs (see Section 3.1).

In all astrophysical runs, we impose a uniform lower magnetic field level of B<sub>0</sub> = 10⁻²⁰ G comoving at z = 60. This extremely low magnetization prevents the formation of spurious numerical effects at the boundary between magnetized and unmagnetized regions in the simulation (in contrast to the primordial models, where there is a non-zero magnetic field everywhere). An overview of the models is given in Table 1.

## 2.3 UHECR simulations

The resulting CMF models used in CRPROPA have a volume of (250 Mpc h⁻¹)<fnref n="3" />, discretized by 1024³ cells of (244 kpc h⁻¹)<fnref n="3" /> volume that contain a uniform field. These models are used to simulate the propagation of UHECRs in the local Universe in order to search for different signatures in the UHECR arrival directions. This is done with CRPROPA 3.0² (Armengaud et al. 2007; Kampert et al. 2013; Batista et al. 2016), a publicly available code to study the propagation of UHECRs. CRPROPA computes all the relevant processes of propagation, this includes Lorentz deflection, energy loss by production of particles and cosmic expansion, photo disintegration and nuclear decay. The code further allows us to track the trajectories of particles in a 3D volume.

We let CRPROPA inject 10⁸ protons with random momentum from random positions. The initial energies range from 1 to 10³ EeV, following a power spectrum of E⁻¹. This choice does not result in the energy spectrum observed in nature, but was used in order to increase number statistics at the highest energies.<fnref n="3" />

<footnote n="2">² https://crpropa.desy.de</footnote>

<footnote n="3">³ A steeper injection spectrum would result in too low accuracy of the measurement of anisotropy around 100 EeV, as can be seen by Eq. 1.</footnote>
<table><caption>**Table 2.** List of the injection models. First column: name of model; second column: set of sources; third column: box length of the simulated volume; fourth column: number of sources in the simulated volume; fifth column: number density of sources.</caption><thead><tr><th>Mnemonic</th><th>Injection scenario</th><th>L<sub>box</sub><br>(Mpc h<sup>-1</sup>)</th><th>N<sub>sources</sub></th><th>n<sub>sources</sub><br>(Mpc<sup>-3</sup>h<sup>3</sup>)</th></tr></thead><tbody><tr><td>homogeneous</td><td>random positions</td><td>250</td><td>10<sup>8</sup></td><td>6,4</td></tr><tr><td>density</td><td>same as homogeneous with p.d.f. = ρ<sub>gas</sub> / Σρ<sub>gas</sub></td><td>250</td><td>10<sup>8</sup></td><td>6,4</td></tr><tr><td>mass halo</td><td>virial haloes, uniform luminosity</td><td>250</td><td>2672</td><td>1.71 · 10<sup>-4</sup></td></tr></tbody></table>

After injection, the energy loss and trajectories of the particles are calculated. In case a trajectory leaves the volume, it is continued on the opposite side. An event is recorded when a trajectory intersects with the observer. This observer is represented by a sphere of radius 800 kpc in the centre of the simulation, which is the defined position of the observer in a constrained simulation. For a discussion on the role of the finite observer size in CRPROPA simulations, we refer the reader to Hackstein et al. (2016).

After intersection, trajectories continue so they may reach another replica of the observer. Environments with strong magnetic fields can trap particles so they arrive at the same observer again. If the same particle is recorded multiple times at the same observer, we randomly chose one of these events. This choice excludes over-counting of trapped particles and no further weighting is necessary.

In a different set of runs, we repeat the process with $10^7$ iron nuclei, taking care also of nuclear decay and disintegration processes, and follow the trajectories of secondary nuclei.

In order to investigate the influence of the distribution of sources, we tested different source models for UHECRs in all the CMF models listed above. In order to bracket the present uncertainties on the degree of isotropy in the distribution of sources, we analyse the extreme case of a homogeneous model, in which we inject each particle at a random position anywhere in the simulated volume. This mimics the absence of structure in the distribution of sources and shows the impact of source distribution in comparison to the other models.

It is generally assumed that sources of UHECRs are powerful sources located in galaxies. Therefore, we assume that the distribution of sources correlates with the LSS. In the *density* model, particles are injected at random positions with a probability density function identical to the gas density, re-normalized by the total gas density in the volume, $p.d.f. = \rho_{gas}/\Sigma\rho_{gas}$. This model with maximum source density reflects a huge number of transient sources that may be found in all types of galaxies, such as gamma-ray bursts or magnetars.

Finally, the *mass halo* model agrees with the lower bounds on source density (~$10^{-4}$ Mpc$^{-3}$, Abreu et al. 2013), where we take as sources the centres of 2672 virial haloes identified in our simulation, each with the same luminosity of UHECRs. This model mimics the case of very few stationary sources, e.g. radio galaxies or active galactic nuclei (AGNs). The precision of the MHD-simulations did not allow us to resolve these structures individually. An overview of the source models can be found in Table 2.

# 3 RESULTS

## 3.1 Simulated extragalactic magnetic fields

Fig. 1 shows the maps of projected gas density (top) and of mean magnetic field along the line of sight (centre, bottom) for the *primordial2R* run (left-hand panels) and for the *astrophysicalR* run (right-hand panels) at $z=0$. Although the different implementa-

tions for gas physics do not significantly change the distribution of gas matter on large scales, the differences in the assumed magnetogenesis scenarios affect the morphological distribution and strength of extragalactic magnetic fields.

In Fig. 2, we present the volume filling factor of the models listed in Table 1. All models have magnetic fields in cluster regions that agree with observational limits. The different *primordial* models show very similar filling factors with dominant strength at ~0.1 nG, close to the upper limit on magnetic field strength in voids from analysis of the CMB anisotropy (Trivedi et al. 2014; Planck Collaboration XIX 2016).

The strong fields in the *astrophysical* models are concentrated in the dense regions of the simulation, which are predominantly filled with very weak fields, at odds with lower limits inferred from the lack of secondary emission around blazar sources (Neronov & Vovk 2010). The filling factors of the *astrophysicalR* and *astrophysical1R* models are almost identical, only in the *astrophysical* model an even smaller volume contains strong fields. Due to the later seeding of magnetic field in all of the astrophysical models, as compared to the primordial models, more of the original, oriented field components survive until $z=0$ and thus a greater influence on the propagation of UHECRs is expected.

## 3.2 Energy spectrum

In Fig. 3, we show the energy spectrum of UHECRs as injected at the sources and measured by the observer. For clarity, the graphs are renormalized by the total number of observed events *N* and multiplied by the inverted energy spectrum at injection, which was set to be $E^{-1}$. Below 100 EeV, the energy spectrum is universal, as predicted by the propagation theorem (Aloisio & Berezinsky 2004). In particular, we find no influence of the underlying magnetic field on the observed energy spectrum, as has been shown in Hackstein et al. (2016).

In the proton injection scenarios, the total number of observed events is *N* = 50 000 with 15 000 events above 10 EeV. In the iron injection scenarios, *N* = 100 000 with 5 000 events above 10 EeV. The fluctuation of these numbers between scenarios with the same initial composition is about 10 per cent. Therefore, number statistics of the observables presented in Sections 3.3 and 3.4 are comparable.

In the proton injection scenarios, the slope above the Greisen-Zatsepin-Kuzmin (GZK) cut-off is not universal but shows significant variation in different source models. We show, both, the Poisson shot-noise as well as the standard deviation for different magnetic field models. They are almost identical, i.e. the error is dominated by statistical fluctuations. Magnetic fields leave no significant impact, as expected for quasi-rectilinear propagation.

The spectrum at ~100 EeV is significantly harder in the *mass halo* injection model, where there is an above-average amount of sources within a few Mpc of the observer. Furthermore, in the *homogeneous* and *mass halo* injection models, protons with up to 800 EeV arrive at the observer in all magnetic field models. However, in the *density*
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Maps of projected gas number density (top) and mean magnetic field along the line of sight (centre, bottom) for the *primordial*2R model (left-hand panel) and for the *astrophysical* R model (right-hand panel) at $z=0$. The gas number density *n* is normalized to the average density in the whole volume, $n/\langle n \rangle$. The magnetic field is shown in µG. Colours are in logarithmic scale. The top and centre panels have a side-length of 200 Mpc $h^{-1}$, the projection axes are the X and Y in the supergalactic coordinates. The bottom panels give a more detailed view on the central 40 Mpc $h^{-1}$. The position of the Milky Way observer considered in this work is exactly at the centre of the box, indicated by a white circle in the top and centre panels. The additional circles show the location of the simulated counterparts of real objects in the local Universe.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Volume filling factor of the models listed in Table 1. The solid lines show the differential filling factor renormalized by 0.1 for clarity, dashed lines show the cumulative filling factor. The grey arrows and shaded area indicate the limits given from observations as listed in the introduction. The yellow line of the *astrophysical*1R model fits exactly with the *astrophysical*1R model.</figcaption></figure>

injection scenario, there are no particles received above 400 EeV. The closer the nearby sources, the higher the number of events that are observed at the most extreme energies and the higher the maximum energy of observed events.

Most particles injected in the simulation never reach the observer and are lost. The injected spectrum plotted in Fig. 3 only shows the injected energies of particles received by the observer. In the proton case, the injected spectrum of observed particles perfectly recreates the injection spectrum used for simulation. In the iron case, multiple secondary nuclei of the same nucleus can reach the observer. In the injected energy spectrum, the primary nucleus is counted once for every secondary nucleus that is observed. This double counting accounts for the sharp increase in the injected spectrum above 40 EeV. At low energies, $\gtrsim 1$ EeV, the injected spectrum is slightly decreased in the stronger magnetic field models. Iron nuclei at low energy are deflected more strongly and are more likely to lose their energy before they reach the observer.

The slope of the observed spectrum is much steeper than in the proton case. The low energies are dominated by the secondary protons of iron injected at the highest energies. Only few events are observed with energies >100 EeV. This is because most of heavy nuclei at those energies disintegrate completely within a few Mpc and distribute their energy evenly among their secondary protons (Epele & Roulet 1998; Allard 2012). Thus, too few events are observed in the iron injection case to measure deviation from isotropy. In conclusion, a sharp cut-off, as observed by extensive air shower arrays (Ivanov 2010; Letessier-Selvon 2014), would hint at a low number of nearby sources or a maximum acceleration energy of protons at the sources that is below the cut-off.

### 3.3 Angular power spectrum

To compute the angular power $C_l$ presented in this section, we first produce full-sky maps of the arrival directions of UHECR events for different minimum energies of considered particles.<fnref n="4" /> These maps are then decomposed into spherical harmonics $\Phi(n) = \sum a_{lm} Y_{lm}(n)$ and $C_l$ is calculated from the ob-

tained amplitudes, $C_l = (2l + 1)^{-1} \sum |a_{lm}|^2$ (cf. Tinyakov & Urban 2015). Finally, the whole spectrum is normalized by the monopole moment, which is $4\pi$ times the square of the average flux.

We present the dipole and quadrupole moment of the angular power spectrum $C_l$ of UHECR arrival directions. These moments were shown to be most promising in the search for anisotropy signals (di Matteo & Tinyakov 2017), but the general trends reported in this section also apply to the octopole moment.

The isotropic prediction is obtained analytically for an isotropic full-sky with $N$ events (Campbell 2015). The mean value of the angular power

$$C_l = 4\pi/N \qquad (1)$$

and the general sample deviation

$$\sigma = \sqrt{\frac{2}{2l + 1} C_l}, \qquad (2)$$

which shows the variation for realizations of a Gaussian random process. For an isotropic sky, both, $C_l$ and $\sigma$, scale with $1/N$. The logarithmic deviation stays constant. In order to account for fluctuation in $C_l$, we show $\sigma$ as error bars for every graph.

Since the value and fluctuation of $C_l$ in an isotropic sky of finite counts are determined by the number of events and our simulations do not reproduce the spectrum observed in nature (cf. Section 3.2), we need to compare to predictions for the simulated spectrum that depends on the injected composition. We indicate with shaded regions the confidence level of anisotropy (C. L. anisotropy). This is obtained from the isotropic prediction and 1, 2, and $3\sigma$ sample deviation, equations (1) and (2). The number of particles $N$ used to calculate the isotropic prediction is the average $N$ observed in each energy bin. The fluctuation of $N$ is about 10 per cent for same injection composition, so the C. L. anisotropy is roughly the same for all models.

In addition to the simulations with magnetic field in all plots, we also present a simulation where the magnetic field is globally set to zero, $B=0$ (black line).<fnref n="5" /> This simulation is shown in order to unambiguously determine the cases where the magnetic field model is important.

We further show the prediction given by a baseline homogeneous model (thick grey line). It shows the average and $1\sigma$ standard deviation of a test group of 27 realizations of a scenario with homogeneous injection in the absence of magnetic fields, $B = 0$. These fully homogeneous scenarios produce the most isotropic results possible in our simulation. The result is not fully isotropic, since it entails all artefacts intrinsic in the simulation, e.g. finite observer effect, overcount of secondary nuclei, and assumed periodicity of the magnetic field and sources (for a detailed discussion see Armengaud, Sigl & Miniati 2005; Hackstein et al. 2016). This makes the homogeneous baseline model a suitable test to find the qualitative contribution of sources and magnetic fields.

The proton injection scenarios are shown in Fig. 4. The prediction from the homogeneous baseline model obtained by the procedure explained above is almost identical to the isotropic prediction. At energies below the GZK-limit of ~40 EeV, the quadrupole angular power is in good agreement with the isotropic prediction for the homogeneous and density injection models presented in the top two panels. In the mass halo injection model, the angular power is above 95 per cent C. L. anisotropy at all energies in virtually all

<footnote n="4">⁴ Due to the hard injection spectrum used in our simulations, the full-sky maps contain too many events at high energies. However, since the observed spectra, in general, are steeper than $E^{-1}$, this effect is negligible.</footnote>

<footnote n="5">⁵ except for the homogeneous plot, where it is given by the homogeneous prediction.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Energy spectrum of UHECRs as injected at the sources (dashed lines) and measured by the observer for a pure proton and a pure iron injected composition (left- and right-hand panels, respectively). The colours and line styles indicate the injection models listed in Table 2. The graphs show the average over all magnetic field models, the standard deviation is indicated by the narrow error bars. The big cyan error bars show the Poisson noise at each second data point. For clarity, the graphs are multiplied by the inverted energy spectrum at injection *E*, and renormalized with the total number of particles *N*.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Angular power $C_l$ of the quadrupole $l = 2$ for all models listed in Table 1 in a pure proton injection scenario. The errorbars indicate sample deviation given by equation (2). From the top to bottom, the panels show the cases of *homogeneous*, *density*, and *mass halo* injection listed in Table 2. The thick grey line is the average and $1\sigma$ standard deviation of the baseline homogeneous model. The shaded regions indicate the 68 per cent, 95 per cent, and 99 per cent C. L. of anisotropy.</figcaption></figure>

of the models. This is in agreement with results from di Matteo & Tinyakov (2017) and Abreu et al. (2013) that show that UHECRs cannot predominantly be protons from few sources in the LSS and that an anisotropic signal should have already been measured for source densities $\lesssim 10^{-4}$ Mpc$^{-3}$.

The magnetic field models do not significantly change the angular power spectrum of arrival directions of UHECR protons at all energies. At very high energies, ~100 EeV, the variation in the coefficients, $C_l$, between the magnetic field models is the lowest, though the number of protons and thus the accuracy is the lowest. The *density* and *mass halo* injection models show a strong deviation from isotropy, whereas the *homogeneous* injection is in good agreement with the prediction from isotropy. The error bars indicate that this feature is not an effect of sample variance, but is statistically significant. This shows that the distribution of nearby sources imposes on the observer an anisotropic signal of UHECRs right below the energy cut-off, where propagation of UHECRs is believed to be quasi-rectilinear. This anisotropic signal can be used to identify the sources of UHECRs.

In the iron injection scenario shown in Fig. 5, almost all models have significantly higher values of $C_l$ below 20 EeV than expected in an isotropic distribution. This energy coincides with $E_{max}/A_{Fe}$, the maximum energy of injected particles $E_{max} = 1000$ EeV divided by the mass number of iron $A_{Fe} = 56$. The predictions from the homogeneous baseline model and the $B=0$ model generally show the highest values. Anisotropy occurs, independent of the source model, due to complete disintegration of heavy nuclei over very short length scales after they have been injected nearby at the highest energies. Due to the high Lorentz-factor, in the absence of deflection, the arrival directions of these secondary nuclei are almost identical, causing an excess of events in direction of the most nearby injection positions (cf. e.g. Lemoine & Waxman 2009). We see that the stronger *primordial* models, in general, show lower $C_l$ values than the weaker *astrophysical* models. The anisotropy produced by the procedure explained above is lowered by CMFs. Since the anisotropy is predominantly produced by nearby sources, only the local field (up to $10 \times$ distance to closest source, Dundović & Sigl 2017) is responsible for this effect. This is in agreement with Sigl et al. (2004), who infer that strong magnetic fields around the observer can suppress large-scale anisotropy.

At the highest energies, ~100 EeV, the number of observed events is too low in the iron injection case to measure the deviation from isotropy.
<figure><img src="image_6.png" /></figure>

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

<figure><img src="image_8.png" /><figcaption>Figure 5. Same as Fig. 4, in a pure iron injection scenario.</figcaption></figure>

During the review process of this manuscript, the Pierre Auger Collaboration reported a significant dipole in the arrival directions of UHECRs with energies >8 EeV at a $5.2\sigma$ level of significance with an amplitude of 6.5 per cent (The Pierre Auger Collaboration et al. 2017) or $C_1 = 0.0050 \pm 0.0025$ in terms of the angular dipole power (Aab et al. 2017). In Figs 6 and 7, we present the dipole moment $l=1$ of the angular power $C_l$ in our simulations and also indicate the recent observation. The features in these graphs are basically the same as discussed for the quadrupole.

Note that the number of particles above 8 EeV in our simulations is different from the amount of particles considered in The Pierre Auger Collaboration et al. (2017). While result is calculated for $\ge 32.000$ events >8 EeV, our simulations have only about $\sim 17.000$ and $\sim 6.000$ events in the proton and iron runs, respectively. Hence, results of our simulations are of lower statistical significance. However, the energy spectra in our simulations are much harder than observed by the Pierre Auger Collaboration, and therefore anisotropic signal from source distribution are expected to be more dominant.

None of the models explored in this paper can reproduce the signal observed in nature with pure proton injection (Fig. 6). Only for

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

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

<figure><img src="image_11.png" /><figcaption>Figure 6. Angular power $C_l$ of the dipole $l=1$ for all models listed in Table 1 in a pure proton injection scenario. The errorbars indicate sample deviation given by equation (2). From the top to bottom, the panels show the cases of homogeneous, density, and mass halo injection listed in Table 2. The thick grey line is the average and $1\sigma$ standard deviation of the baseline homogeneous model. The shaded regions indicate the 68 per cent, 95 per cent, and 99 per cent C. L. of anisotropy. The red point corresponds to the amplitude of the recent dipole signal reported by Auger.</figcaption></figure>

the mass halo model there is a small overlap of $1\sigma$ deviations with the Auger measurement. The level of anisotropy does not decrease strongly from 8 to 4 EeV. This indicates that a strong dipole in the distribution of nearby sources is necessary to reproduce the Auger signal with a light injection composition. In that scenario, the dipole angular power $C_1$ increases at most linear with energy between 1 and 10 EeV in our simulations. The amplitude is proportional to the square root of $C_1$ and increases too slow compared to observations by the Pierre Auger Observatory. This makes a light injection composition of UHECRs at the highest energies unlikely in view of the recent observation.

Injection of iron nuclei, Fig. 7, results in a dipole similar to that observed by Auger – in amplitude, not in significance. Further, $C_1$ increases roughly quadratically for heavy injection composition
<figure><img src="image_12.png" /></figure>

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

<figure><img src="image_14.png" /><figcaption>**Figure 7.** Same as Fig. 6, in a pure iron injection scenario.</figcaption></figure>

and compares better to results of Pierre Auger Collaboration. In the *homogeneous* model, magnetic fields can suppress the signal to agree with the isotropic prediction. In the *density* and *mass halo* model, the anisotropic signal is dominated by the distribution of sources and not suppressed efficiently by magnetic fields.

Our results suggest that the dipole signal in UHECRs observed by the Pierre Auger Observatory may be the product of clustering of secondary nuclei in direction of the nearby sources.

## 3.4 Composition

In Fig. 8, we show the average mass number $\langle A \rangle$ of observed events as function of energy in the iron injection scenarios. At low energies, the composition is very light since secondary protons of injected iron nuclei dominate observations. At $\lesssim 20$ EeV $\approx E_{max}/Z_{Fe}$, there is a steep increase in $\langle A \rangle$. This coincides with the maximum energy of secondary protons. At higher energies, only the (partly disintegrated) primary nuclei are observed. All magnetic field and injection models show a very similar slope of $\langle A \rangle$. We conclude that CMFs in agreement with observational upper limits, in general, are too weak

<figure><img src="image_15.png" /><figcaption>**Figure 8.** Average mass number $\langle A \rangle$ of UHECRs observed at different energies. The errorbars show the $1\sigma$ standard deviation. The colours indicate the magnetic field model listed in Table 1 and the linestyle shows the source model listed in Table 2.</figcaption></figure>

to impose a significant difference in the all-sky average composition of UHECRs.

## 4 CONCLUSIONS

We have studied the influence of CMFs on the propagation of UHECRs using MHD-simulations with different models for seeding of magnetic fields for both, primordial and astrophysical, processes. We found no evidence that magnetic field seeding scenarios could be distinguished via the use of the angular power spectrum of the spherical harmonics decomposition of the full-sky of arrival directions of UHECRs.

We have studied the influence of different source scenarios on the energy spectrum of UHECRs and on the angular power of anisotropy. We have found that for a pure proton composition, the slope of the energy spectrum at energies > 100 EeV depends on the number of, and distance to, the most nearby sources. The closer the sources, the harder the energy spectrum. If only iron is injected, almost no events are observed above that energy. Thus, the sharp energy cut-off observed with extensive air shower arrays (Ivanov 2010; Letessier-Selvon 2014) might suggest a low number of sources in the near vicinity of the observer if the cut-off does not coincide with the maximum energy of proton acceleration.

We have investigated the angular power spectrum of arrival directions. We have found that there is a clear deviation from isotropy, $\lesssim 100$ EeV, if the distribution of sources follows the LSS. This offers the chance to identify the sources with future full-sky measurements (Dawson, Fukushima & Sokolsky 2017) and high number statistics at the highest energies.

We were able to reproduce the dipole in the arrival directions of UHECRs > 8 EeV recently reported by the Pierre Auger Collaboration (The Pierre Auger Collaboration et al. 2017) with all our source models, but only using pure iron injection composition instead of protons. Our results indicate that the observed dipole is the result of clustering in direction of nearby sources of heavy nuclei (Lemoine & Waxman 2009). Strong magnetic fields might be necessary to explain the absence of anisotropy signal in the higher multipoles. Exploring such possibilities (also joined with a more thorough exploration of the role of UHECRs composition in the production of a dipole excess) will be subject of forthcoming work.

For the injection of protons from the virial haloes with a very low number density, around the limit from Abreu et al. (2013, $\sim 10^{-4}$ Mpc$^{-3}$), we have found 95 per cent C. L. quadrupolar anisotropy at all energies, in conflict with present observations (Pierre Auger Collaboration 2012; Aab et al. 2014). This confirms
the findings of di Matteo & Tinyakov (2017) that UHECRs cannot primarily be protons from few sources in the LSS.

Finally, we have analysed the observed composition of UHECRs via the average mass number of events. There is no evidence that CMFs significantly influence the all-sky composition of UHECRs at all energies.

In our study, we did not account for the influence of the magnetic field of the Milky Way, but energy losses are negligible on galactic scales. Furthermore, the angular power spectrum at large scales has been shown to have low impact of deflections in the Galactic magnetic field (Tinyakov & Urban 2015; di Matteo & Tinyakov 2017).

In summary, with newer constrained simulations of the local Universe, we confirmed our previous findings (Hackstein et al. 2016), i.e. that the properties of observed UHECRs do not seem to carry much information on the genesis and distribution of extragalactic magnetic fields. This in turn strengthens the possibility of performing 'UHECRs astronomy' (Dolag et al. 2004), thus motivating further investigations on the origin of UHECRs across a wide range of energies where the impact of the Galactic magnetic field should be sub-dominant.

## ACKNOWLEDGEMENTS

Computations described in this work were performed using the ENZO code (http://enzo-project.org), which is the product of a collaborative effort of scientists at many universities and national laboratories. The constrained initial conditions have been developed in the context of the Constrained Local UniversE Simulations (CLUES) project (https://www.clues-project.org/). FV acknowledges financial support from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska-Curie grant agreement no. 664931 and the ERC Starting Grant MAGCOW, no. 714196. JGS acknowledges support from the Astronomy ESFRI and Research Infrastructure Cluster ASTERICS project, funded by the European Commission under the Horizon 2020 Programme (GA 653477). We acknowledge the usage of computational resources on the JURECA cluster at the Juelich Supercomputing Centre (JSC), under projects no. 10755, 11823, 9016 and 8998, and on the Piz-Daint supercluster at CSCS-ETHZ (Lugano, Switzerland) under project s701. SH, FV, and MB additionally thank G. Sigl and A. Dundovic for fruitful discussion as well as the other developers of CRPROPA for making the code public. We further acknowledge the use of computational resources at the Rechenzentrum of the University of Hamburg. The whole group would also like to thank the referee for fruitful comments and careful reading of the manuscript. We also thank S. Mollerach for adjuvent comments and discussion on the manuscript.

## REFERENCES

* Aab A. et al., 2014, ApJ, 794, 172

* Aab A. et al., 2017, J. Cosmol. Astropart. Phys., 6, 026

* Abreu P. et al., 2013, J. Cosmol. Astropart. Phys., 1305, 009

* Allard D., 2012, Astropart. Phys., 39, 33

* Aloisio R., Berezinsky V., 2004, ApJ, 612, 900

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

* Araya-Melo P. A., Aragón-Calvo M. A., Brüggen M., Hoeft M., 2012, MNRAS, 423, 2325

* Armengaud E., Sigl G., Miniati F., 2005, Phys. Rev., D, 72, 043009

* Armengaud E., Sigl G., Beau T., Miniati F., 2007, Astropart. Phys., 28, 463

* Batista R. A. et al., 2016, J. Cosmol. Astropart. Phys., 2016, 038

* Beck R., 2016, A&AR, 24, 4

* Bernet M. L., Miniati F., Lilly S. J., 2013, ApJ, 772, L28

* Bertschinger E., 1987, ApJ, 323, L103

* Bistolas V., Hoffman Y., 1998, ApJ, 492, 439

* Blasi P., 2013, A&AR, 21, 70

* Bonafede A., Vazza F., Brüggen M., Murgia M., Govoni F., Feretti L., Giovannini G., Ogrean G., 2013, MNRAS, 433, 3208

* Broderick A. E., Chang P., Pfrommer C., 2012, ApJ, 752, 22

* Brown S. D., 2011, J. Astrophys. Astron., 32, 577

* Brown S. et al., 2017, MNRAS, 468, 4246

* Brüggen M., Ruszkowski M., Simionescu A., Hoeft M., Dalla Vecchia C., 2005, ApJ, 631, L21

* Bryan G. L. et al., 2014, ApJS, 211, 19

* Campbell S. S., 2015, MNRAS, 448, 2854

* Das S., Kang H., Ryu D., Cho J., 2008, ApJ, 682, 29

* Dawson B. R., Fukushima M., Sokolsky P., 2017, preprint (arXiv:1703.07897)

* Dedner A., Kemm F., Kröner D., Munz C.-D., Schnitzer T., Wesenberg M., 2002, J. Comput. Phys., 175, 645

* di Matteo A., Tinyakov P., 2017, MNRAS, preprint (arXiv:1706.02534)

* Dolag K., 2006, Astron. Nachr., 327, 575

* Dolag K. M., Bartelmann M., Lesch H., 1999, A&A, 348, 351

* Dolag K., Grasso D., Springel V., Tkachev I., 2004, J. Exp. Theor. Phys. Lett., 79, 583

* Donnert J., Dolag K., Lesch H., Muller E., 2009, MNRAS, 392, 1008

* Dova M. T., 2016, preprint (arXiv:1604.07584)

* Dundović A., Sigl G., 2017, preprint (arXiv:1710.05517)

* Epele L. N., Roulet E., 1998, J. High Energy Phys., 10, 009

* Feretti L., Giovannini G., Govoni F., Murgia M., 2012, A&AR, 20, 54

* Ganon G., Hoffman Y., 1993, ApJ, 415, L5

* Hackstein S., Vazza F., Brüggen M., Sigl G., Dundovic A., 2016, MNRAS, 462, 3660

* Harari D., Mollerach S., Roulet E., 2000, J. High Energy Phys., 02, 035

* Harari D., Mollerach S., Roulet E., 2002a, J. High Energy Phys., 7, 006

* Harari D., Mollerach S., Roulet E., Sánchez F., 2002b, J. High Energy Phys., 3, 045

* Heß S., Kitaura F.-S., Gottlöber S., 2013, MNRAS, 435, 2065

* Hillas A. M., 1984, Annu. Rev. Astron. Astrophys., 22, 425

* Hoffman Y., Ribak E., 1991, ApJ, 380, L5

* Hoffman Y., Ribak E., 1992, ApJ, 384, 448

* Ivanov A. A., 2010, ApJ, 712, 746

* Jasche J., Wandelt B. D., 2013, MNRAS, 432, 894

* Kampert K.-H., Kulbartz J., Maccione L., Nierstenhoefer N., Schiffer P., Sigl G., van Vliet A. R., 2013, Astropart. Phys., 42, 41

* Kim K. S., Lilly S. J., Miniati F., Bernet M. L., Beck R., O'Sullivan S. P., Gaensler B. M., 2016, ApJ, 829, 133

* Kitaura F.-S., 2013, MNRAS, 429, L84

* Klypin A., Hoffman Y., Kravtsov A. V., Gottlöber S., 2003, ApJ, 596, 19

* Kotera K., Lemoine M., 2008, Phys. Rev., D, 77, 023005

* Kravtsov A. V., Klypin A., Hoffman Y., 2002, ApJ, 571, 563

* Kurganov A., Tadmor E., 2000, J. Comput. Phys., 160, 241

* Lavaux G., 2010, MNRAS, 406, 1007

* Lavaux G., Mohayaee R., Colombi S., Tully R. B., Bernardeau F., Silk J., 2008, MNRAS, 383, 1292

* Lemoine M., Waxman E., 2009, J. Cosmol. Astropart. Phys., 11, 009

* Letessier-Selvon A., 2014, Braz. J. Phys., 44, 560

* Neronov A., Vovk I., 2010, Science, 328, 73

* Pierre Auger Collaboration, 2012, ApJS, 203, 34

* Pshirkov M. S., Tinyakov P. G., Urban F. R., 2016, Phys. Rev. Lett., 116, 191302

* Planck Collaboration XVI, 2014, A&A, 571, A16

* Planck Collaboration XIX, 2016, A&A, 594, A19

* Ryu D., Schleicher D. R. G., Treumann R. A., Tsagas C. G., Widrow L. M., 2012, Space Sci. Rev., 166, 1

* Shu C.-W., Osher S., 1988, J. Comput. Phys., 77, 439

* Sigl G., Miniati F., Ensslin T. A., 2003, Phys. Rev., D68, 043002

* Sigl G., Miniati F., Ensslin T., 2004, Nucl. Phys. Proc. Suppl., 136, 224

* Sigl G., Miniati F., Enßlin T. A., 2004, Phys. Rev. D, 70, 043007
Sorce J. G., 2015, MNRAS, 450, 2644

Sorce J. G., Courtois H. M., Gottlöber S., Hoffman Y., Tully R. B., 2014, MNRAS, 437, 3586

Sorce J. G. et al., 2016, MNRAS, 455, 2078

Stanev T., 1997, AJ, 479, 290

Takami H., Sato K., 2008, Astrophys. J., 681, 1279

Takami H., Inoue S., Yamamoto T., 2012, Astropart. Phys., 35, 767

The Pierre Auger Collaboration et al., 2017, preprint (arXiv:1709.07321)

Tinyakov P. G., Tkachev I. I., 2005, Astropart. Phys., 24, 32

Tinyakov P. G., Urban F. R., 2015, Sov. J.Exp. Theor. Phys., 120, 533

Trivedi P., Subramanian K., Seshadri T. R., 2014, Phys. Rev. D, 89, 043523

Tully R. B. et al., 2013, AJ, 146, 86

Vallée J. P., 2004, New Astron. Rev., 48, 763

van de Weygaert R., Bertschinger E., 1996, MNRAS, 281, 84

van Vliet A. R., 2014, PhD thesis, Hamburg University

Vazza F., Dolag K., Ryu D., Brunetti G., Gheller C., Kang H., Pfrommer C., 2011, MNRAS, 418, 960

Vazza F., Ferrari C., Brüggen M., Bonafede A., Gheller C., Wang P., 2015, A&A, 580, A119

Vernstrom T., Gaensler B. M., Brown S., Lenc E., Norris R. P., 2017, MNRAS, 467, 4914

Wang P., Abel T., Kaehler R., 2010, New Astron., 15, 581

Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94

Yoshiguchi H., Nagataki S., Tsubaki S., Sato K., 2003, AJ, 586, 1211, (Erratum: Yoshiguchi H., Nagataki S., Tsubaki S., Sato K., 2004, AJ, 601, 592)

## APPENDIX A: RE-WEIGHTED ENERGY SPECTRUM

The energy spectra presented in this work do not recreate the spectrum observed in nature. This can be achieved by modifying the injection spectrum, in particular, by using a softer spectral index and

introducing an exponential cut-off. The resulting injection spectrum has the form

$$ \frac{dN}{dE} \propto E_0^{-\gamma} e^{-E_0/(Z_0 R_{\text{max}})}, \qquad (A1) $$

with initial energy $E_0$, initial charge number $Z_0$, spectral index $\gamma$, and maximum rigidity $R_{\text{max}} = E_{0,\text{max}}/Z_0$. The modification can be done in post-processing by multiplying every event with a specific weight factor (Armengaud et al. 2005; van Vliet 2014)

$$ w(E_0, Z_0) = E_0^{\gamma_{\text{init}}-\gamma} e^{-E_0/(Z_0 R_{\text{max}})}, \qquad (A2) $$

where $\gamma_{\text{init}}$ is the spectral index used for the simulation.

In order to obtain the correct spectrum at injection, we fit the observed spectral index between the ankle $E_A \approx 5$ eV and the cut-off $E_C \approx 20$ eV, which is observed to be $\gamma = 2.63 \pm 0.04$ (Letessier-Selvon 2014). The best fit in the proton injection scenarios is an injection index of $\gamma = 2$, as expected for Fermi acceleration. The best fit for the iron injection scenarios is $\gamma = 2.4$.

The shape of the spectrum beyond the cut-off energy $E_C$ is recreated well by using $E_{0,\text{max}} = 100$ eV for the exponential cut-off. The maximum rigidity is then $R_{\text{max,p}} = 100$ EV for the proton injection scenarios and $R_{\text{max,Fe}} = 100/26$ EV $\approx 3.8$ EV for iron injection scenarios. The resulting spectra are shown in Fig. A1.

After re-weight, the effective number of observed particles is $N_{\text{eff,p}} \approx 7000$ in the proton injection scenarios and $N_{\text{eff,Fe}} \approx 600$ in the iron injection scenarios. The isotropic prediction for the angular power spectrum depends on the number of particles (see eq. 1). Therefore, after re-weight the isotropic prediction increases everywhere by about an order of magnitude at least. Accordingly, the colour bands are raised in Figs 4–7. All re-weighted scenarios are below 68 per cent C. L. anisotropy at all energies.

<figure><img src="image_16.png" /><figcaption>**Figure A1.** Re-weighted energy spectrum of UHECRs as injected at the sources (dashed lines) and measured by the observer for a pure proton and a pure iron injected composition (left- and right-hand panels, respectively). The colours and line styles indicate the injection models listed in Table 2. The graphs show the average over all magnetic field models, the standard deviation is indicated by the narrow error bars. The big cyan error bars show the Poisson noise at each second data point. The graphs are renormalized with the total number of particles $N$ and multiplied by $E^{-3}$ to enable better comparison to the figures presented in Letessier-Selvon (2014).</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.