Impact of the Magnetic Horizon on the Interpretation of the Pierre Auger Observatory Spectrum and Composition Data

. The flux of ultra-high energy cosmic rays reaching Earth above the ankle energy (5 EeV) can be described as a mixture of nuclei injected by extragalactic sources with very hard spectra and a low rigidity cutoff. Extragalactic magnetic fields existing between the Earth and the closest sources can affect the observed CR spectrum by reducing the flux of low-rigidity particles reaching Earth. We perform a combined fit of the spectrum and distributions of depth of shower maximum measured with the Pierre Auger Observatory including the effect of this magnetic horizon in the propagation of UHECRs in the intergalactic space. We find that, within a specific range of the various experimental and phenomenological systematics, the magnetic horizon effect can be relevant for turbulent magnetic field strengths in the local neighbourhood in which the closest sources lie of order B rms ≃ (50 − 100) nG (20 Mpc / d s )(100 kpc / L coh ) 1 / 2 , with d s the typical intersource separation and L coh the magnetic field coherence length. When this is the case, the inferred slope of the source spectrum becomes softer and can be closer to the expectations of diffusive shock acceleration, i.e., ∝ E − 2 . An additional cosmic-ray population with higher source density and softer spectra, presumably also extragalactic and dominating the cosmic-ray flux at EeV energies, is also required to reproduce the overall spectrum and composition results for all energies down to 0.6 EeV.

To understand the properties of the sources of cosmic rays (CRs) required in order to account for the spectrum and composition inferred from data collected at the Pierre Auger Observatory, combined fits to the measurements of the CR flux and distribution of the depth of shower maximum (X max ) were performed considering different astrophysical source scenarios [1,2] (see also [3][4][5][6]).These fits adopted continuous distributions of CR sources, eventually allowing for a redshift evolution of their emissivities.The emitted particles were then propagated including the attenuation effects due to the CR interactions with the background radiation, and the resulting fluxes at the Earth were compared to observations in order to obtain the best fitting source parameters.One source population dominates the fluxes above the ankle energy (∼ 5 EeV), while a second source population is required to also explain the observations below the ankle, as in particular shown in [2][3][4], with each population having sources with specific spectral and composition properties.Since the observations indicate that the CR composition becomes increasingly heavier above the ankle energy [7,8], this can naturally result if the different mass components have cutoffs which depend on the rigidity R = E/Z, i.e. cutoff energy proportional to the atomic number Z, as expected from electromagnetic acceleration processes 1 , so that the light components do not reach the highest energies.Moreover, in order for the flux of the heavier components dominating at the highest energies to be sufficiently suppressed below the ankle, where a light composition is inferred, the individual spectral shapes of the different elements contributing to the high-energy population need to be quite hard.In particular, considering that below the cutoffs the source spectra have a power-law shape E −γ , values of γ < 1 are inferred (and in some cases as low as −2) for the high-energy population, which is in tension with the value γ ≃ 2 that is expected from diffusive shock acceleration (DSA).Similar values for the spectral index of the high-energy population were obtained when fitting the spectrum, composition and arrival direction distribution at energies above 10 19 eV [9].The energy dependent magnetic confinement of heavy nuclei and their photodisintegration in the source environment has been considered as a possible explanation of the suppression of the flux of the high-energy component at low energies in [4,[10][11][12][13][14][15][16], leading to an emission spectrum harder than that produced in the acceleration process.The population dominating the flux below the ankle is instead required to have a very soft spectrum (γ > 3) and a mix of protons and intermediate-mass nuclei.The very soft spectrum may eventually be due to the superposition of many sources with a harder spectrum and a distribution of cutoff rigidities [17].The proton component below the ankle could alternatively result from interactions in the source environment of the intermediate and heavy-mass nuclei from the high-energy population, producing nucleons, among which neutrons that can escape the magnetized medium and decay into protons in flight.This model has been tested already in [2,18].The previous combined fits of the spectrum and composition [1,2] have considered a continuous distribution of CRs sources in the universe.In the more realistic case of a discrete source distribution, the presence of sizeable extragalactic magnetic fields between the closest sources and the Earth can affect the shape of the spectrum of cosmic rays reaching us.In particular, an alternative explanation for the hardness of the high-energy population spectrum is that it is due to a magnetic horizon effect which, as a consequence of the CR diffusive propagation through the intergalactic turbulent magnetic fields, could suppress the flux at low energies [19][20][21].This suppression is relevant when the time for CR diffusion from the closest sources becomes longer than the age of the sources, so that the low-energy cosmic rays have not enough time to reach the Earth, and for this to happen one then needs strong magnetic fields as well as relatively large intersource distances (i.e.small source densities).Following the ideas in [19,[21][22][23][24], we here extend the combined fit analysis of the spectrum and composition data by including the magnetic horizon effects associated to the finite intersource separations and the presence of extragalactic magnetic fields, which are two ingredients that should be included in a more realistic description of reality.We explore in particular whether the resulting suppression of the flux at low energies can allow for softer source spectra so as to alleviate the existing tension with the expectations from DSA.We don't consider the Galactic magnetic field since it is not expected to significantly modify the spectrum of extragalactic cosmic rays.

The source populations
The general framework we consider in this study to reproduce the energy spectrum both below and above the ankle feature is that already explored in [2], consisting of two different source populations.One source population dominates the CR flux above a few EeV while a second one dominates at lower energies.For each of these populations the source spectra of the different mass components are considered to be described as power laws having a rigidity dependent cutoff which strongly suppresses the fluxes above an energy ZR cut .This kind of cutoffs are expected when the acceleration has an electromagnetic origin, so that the maximum energy achieved is proportional to the charge of the particle, and when there are no strong interactions in the acceleration region. 2 In that case, the differential particle generation rate of each component of atomic number Z and mass number A, per unit volume, energy and time, is with a = (L, H) identifying the population dominating at low and high energies respectively.For each population, the normalization Qa 0 is the present total differential rate of CR emission per unit energy, volume and time, at the reference energy E 0 (smaller than the hydrogen cutoff R a cut and taken here as 1 EeV), at which the relative source fractions of the different elements are f a A .We consider that five representative elements are emitted at the sources: H, He, N, Si and Fe.
The factor ξ a (z) parameterises the evolution of the emissivity as a function of the redshift z, for which we consider here two possibilities.The first is that of non-evolving sources (NE) having ξ a NE (z) = 1, and the other assumes that the emissivities scale with the star formation rate (SFR) as parameterised in [26], i.e.
with a steep decline at higher redshifts.We will actually consider maximum redshifts z max = 1 for the NE case and z max = 4 for the SFR case, since the contribution from CR sources farther away is negligible.For the cutoff function F cut we adopt the following parametrization [2] F cut (y) = sech(y ∆ ), with the parameter ∆ determining the steepness of the cutoff shape, which turns out to have a significant impact on the fit.We will consider the three representative values ∆ = 1, 2 and 3. 3The particles emitted at the sources are propagated up to the Earth using the SimProp v2r4 software [27].The resulting fluxes depend on the nuclear photo-disintegration cross sections, as well as on the extragalactic background light model considered to evaluate the interactions during propagation.We adopt in the analysis the photodisintegration cross section from TALYS [28] and the extragalactic background radiation from Gilmore et al. [29].The inferred composition depends on the hadronic model used to interpret the X max measurements, which was found in [1,2] to be the assumption that mostly affected the values of the fitted fractions.We will explore the dependence of the results on this by considering both EPOS-LHC [30] and Sibyll 2.3d [31] hadronic interaction models.
3 The magnetic horizon effect An alternative explanation [21] for the hardness of the observed high-energy component spectra is that it is not due to the hardness of the source spectra, but to a magnetic horizon effect [32][33][34] which, as a consequence of the CR diffusive propagation through the intergalactic turbulent magnetic fields, could increasingly suppress their flux for decreasing energies.We will here consider the case of steady sources, and given the large z max values considered the characteristic time for the source emission is of the order of the age of the universe.Note that the suppression effects at low energies could be more pronounced if the sources were transient rather than steady, in which case smaller strength of the magnetic fields and/or smaller intersource distances would be required, but in this case the results will depend on the specific emission histories of the different nearby CR sources and their distances [35,36].
Magnetic fields are known to permeate the Universe, having different strength on different scales, and they hence affect the propagation of the charged cosmic rays.In our galaxy they have a strength of several µG and extend over scales up to tens of kpc.They can be modelled using several components, such as the regular ones in the disk and in the halo which are coherent over kpc scales, or the random and striated turbulent components which have a typical coherence length of about 50 pc (see e.g.[37]).These components will deflect the incoming UHECR and also lead to a diffusive behaviour for the CRs with energies below about 0.1Z EeV.In clusters of galaxies the turbulent magnetic fields can reach values of 1 to 10 µG, with typical coherence lengths of 1 to 100 kpc [38,39].In large-scale structure filaments the magnetic fields are more uncertain, but may well have strengths in the 10 to few 100 nG range, with coherence lengths in the 10 kpc to 1 Mpc range [39][40][41], while in the voids of the large-scale structures they are expected to be much weaker, with strengths smaller than 1 nG [42,43].The strength of these fields depends on the mechanism producing them.In particular, they could result from the amplification via flux conserving gravitational compression of primordial seeds, such as those left over from an inflationary period or those produced in phase transitions in the early universe, or alternatively they could result from outflows of galactic fields that were amplified by some dynamo process.
The presence of extragalactic magnetic fields (EGMF) can be the dominant contributor to a magnetic horizon effect and may hence affect the observed cosmic-ray spectrum.In particular, EGMF are expected to be enhanced in our local neighbourhood where the closest UHECR sources should lie, such as within the Local Supercluster which has a typical radial extent of about 20 Mpc.Also note that the assumption of equipartition between the energy in thermal gas motion and in the magnetic field leads to an estimation of the magnetic field strength within filaments and sheets of the large scale structure of order 10 2 nG [44].We will model for simplicity the EGMF as being turbulent and isotropic, parameterized by its root-mean-square amplitude (B rms ) and coherence length (L coh ), considering a Kolmogorov spectrum for the turbulence.A critical energy can be defined as that for which the effective Larmor radius associated to B rms equals the coherence length, and for a particle of atomic number Z it is given by E crit ≡ |e|ZB rms L coh ≡ ZR crit , with the critical rigidity being the critical energy for a H nucleus.The critical energy separates two different propagation regimes.The first one is that of resonant diffusion, at energies E ≪ E crit , in which deflections are large even before traversing a distance L coh .The second one is for energies E ≫ E crit , in which case the propagation is instead quasirectilinear within a coherence length and the diffusive regime is only attained after traversing much longer distances.
We will consider a distribution of UHECR sources with uniform density n s , with the average distance between them being d s = n −1/3 s .When the particle rigidities are low enough so that the CR diffusive travel time from the closest sources in the presence of an EGMF becomes larger than the source age, their flux will get significantly suppressed.It is useful to quantify this effect through a suppression function G(E) ≡ J(E)/J(E) ds→0 , given by the ratio between the actual flux at Earth from the discrete source distribution to the flux that would result in the limit of a continuous source distribution (with the same emissivity per unit volume). 4Using simulations of the propagation of particles in turbulent magnetic fields with a Kolmogorov spectrum performed in an extended implementation of the SimProp code where both interactions with radiation fields and magnetic deflections are accounted for, the suppression function has been parameterised as [45] where x ≡ E/(ZR crit ) and X s = d s / √ r H L coh is the normalized intersource distance, with r H = c/H 0 the Hubble radius (in terms of the speed of light c and the present day Hubble constant H 0 ≃ 70 km s −1 Mpc −1 ), so that The parameters a, b, α and β are sensitive to the distribution of the initial redshifts of the particles that reach the Earth.They hence depend on the assumed cosmological evolution of the source population emissivity as well as on whether the particles are primaries emitted at the sources or secondaries produced by photo-disintegration interactions during their propagation.The suppression also depends slightly on the spectral index of the sources, and the values of the parameters obtained in [45] are tabulated in the Appendix A, both for the NE and SFR scenarios.
The spectra of the different mass components reaching the Earth in the presence of EGMF can be obtained as the product of those in the absence of magnetic fields times the corresponding suppression factor G. Thus, the magnetic horizon effect is accounted through two parameters: the critical rigidity R crit and the normalized intersource distance X s .This description allows us to probe a wide range of values of the magnetic field amplitudes and coherent lengths as well as of the source densities.It is important to keep in mind that due to the diffusion the contribution from the faraway sources becomes strongly suppressed for decreasing rigidities and eventually the flux at low energies is dominated by that from the nearby sources.Thus, the magnetic field which is relevant for the suppression due to the magnetic horizon effect is the one between the closest sources and the observer, with the closest source being on average at a distance of 0.55d s [21].If one considers magnetic field amplitudes in the range 1 nG ≲ B rms ≲ 200 nG and coherence length such that 25 kpc ≲ L coh ≲ 1 Mpc, one expects that 0.022 EeV ≲ R crit ≲ 180 EeV.Moreover, if the intersource average distance is in the range 4 Mpc ≲ d s ≲ 40 Mpc one should have that 0.05 ≲ X s ≲ 4. We will hence consider parameters R crit and X s within these ranges when performing the fits that include the magnetic horizon effects.
The magnetic suppression factor G is displayed in Fig. 1 as a function of E/E crit , for different values of X s and for the NE (left) and SFR (right) evolutions, considering γ = 1.The implied strong hardening of the spectrum at low energies is apparent, with the suppression occurring in general in the regime of resonant diffusion, i.e. for energies below E crit ≡ ZR crit .This effect is clearly rigidity dependent, given its magnetic origin.One also finds that for increasing X s the suppression appears for larger values of E/E crit .Hence, increasing X s and simultaneously decreasing E crit the different curves can lead to comparable suppressions as a function of the energy, what will lead to some approximate degeneracies between these parameters, although the shape of the suppression becomes steeper for increasing values of X s .For the same values of the parameters X s and E crit the suppressions are milder for the SFR case than for the NE one, given that due to the enhanced emission at high redshifts the particles have on average more time to arrive from their sources in the first case (note that for diffusing particles the redshift is a measure of the time travelled rather than of the source distance).We show the results both for primaries (solid lines), secondary protons (dashed lines) and intermediate secondary nuclei (dotted lines), where the primaries are those in which the detected nucleus is in the same mass group as the emitted one, considering the different mass groups as those with values of A of 1 (H), 2-4 (He), 5-16 (N), 17-30 (Si) and 31-56 (Fe).The intermediate secondary nuclei are those in a lighter mass group than the primary injected particle.The milder suppression of the secondaries is also understood because their average redshift of production in photodisintegration processes is higher than the average production redshift of the primaries, given that the photodisintegrations processes get enhanced at higher redshifts.Also note that the suppression of the secondary protons is quite similar to the suppression of the secondary nuclei of the same rigidity.
In the scenarios with two source populations that we are considering, we will assume that the low-energy component providing the bulk of the flux at few EeV energies arises from sources which are more abundant than those of the high-energy component, having then a smaller intersource distance (smaller X s ).The magnetic suppression will in this case affect the LE component at energies much lower than the HE component [32,46], and we will hence neglect the magnetic horizon effects on the low-energy component.Since this approximation may become less accurate for decreasing energies and, in addition, a non-negligible Galactic contribution of heavy elements may also extend above the second-knee feature present at 0.1 EeV, the model predictions at energies below 1 EeV are expected to be affected by those effects.
4 Combined fit to the spectrum and composition

Flux model
The flux of particles reaching the Earth from each population of sources a can then be obtained for a given model of the sources and extragalactic magnetic field parameters as with dt dz = 1 where E ′ and E denote the energy observed at Earth and that at the sources, respectively, Ω Λ ≃ 0.7 and Ω m ≃ 0.3.The effect of the interactions with background photons is accounted for by dη A ′ ,A (E ′ , E, z)/dE ′ , which represents the differential probability that particle with energy E ′ and mass number A ′ reaches z = 0 when a particle with energy E and mass A is injected at a redshift z, and this quantity is obtained from simulations computed using SimProp.

Fit procedure
The parameters of the model are obtained by maximising the likelihood so that the assumed model reproduces the observed data.To do so, we follow the steps outlined in [1,2].The test statistic used to parameterise the goodness of fit is the deviance, which is defined as twice the negative logarithm of the likelihood ratio between each model and one that would perfectly describe the observed data (called the saturated model ), D = −2 ln (L/L sat ).This implies that for the parameters at which the deviance is at a minimum the likelihood of the model is maximized.Since both the energy spectrum and X max distributions were measured independently, the complete likelihood of our model is simply the product L J • L Xmax .Here L J represents the likelihood of the energy spectrum, defined as where i denotes the i-th detection energy bin; J obs i and J mod i are the observed flux and the one predicted by the model, respectively, and σ i the associated uncertainties in the measurements.Meanwhile, L Xmax represents the likelihood of the X max distribution measurements, following the expression where k obs i,j is the number of events in the i-th energy bin and the j-th bin of the X max distribution, n obs i is the total number of event in the i-th energy bin, and G mod i,j are the corresponding model predictions obtained from the modified Gumbel functions which account for the detection and resolution effects as described in [2].
The deviance is minimized using the Minuit library [47].The uncertainties in the spectral parameters, fractions and magnetic horizon parameters are obtained via the MINOS procedure included in the Minuit package.These correspond to the change in each parameter for which the deviance increases by one unit when minimising with respect to the rest of the parameters.The uncertainties in the source emissivities are obtained via Monte Carlo simulations in which the different parameters are allowed to vary within their uncertainties.Table 1.Parameters of the fit to the spectrum and Xmax distributions in the absence of a magnetic horizon effect, for the EPOS-LHC and Sibyll 2.3d hadronic interaction models and no source evolution, with the corresponding deviance D and number of fitted data points N .Cutoff shapes with ∆ = 1, 2 and 3 are considered.

Data sets
We fit the energy spectrum determined by the Pierre Auger Observatory using the events detected by the Surface Detector array.The array with stations separated by 1500 m is used above 2.5 EeV while the denser array with stations separated by 750 m is used for smaller energies [48].The energy range considered covers from 10 17.8 eV to 10 20.2 eV, in logarithmic bins of width ∆ log 10 (E/eV) = 0.1.For the composition we fit the X max distributions measured using the Fluorescence Detector telescopes in the energy range from 10 17.8 eV to 10 19.6 eV in logarithmic bins of width ∆ log 10 (E/eV) = 0.1, plus one additional bin including the events with energies above 10 19.6 eV [7].The X max distributions are binned in bins of width ∆X max = 20 g cm −2 in each energy interval.These are the same data sets considered in [2,24].

Results obtained in the absence of magnetic fields
As a starting point, we first present in Table 1 the results of the fit performed in the absence of magnetic fields and for non-evolving sources, for the different cutoff functions considered and for both EPOS-LHC and Sibyll 2.3d hadronic interaction models, in line with the analysis in [2], and with compatible results. 6 The spectrum parameters (γ and R cut ) of the high-energy and low-energy populations are reported, as well as the obtained deviance.It is interesting to note the effect of the cutoff function on the spectral parameters of the HE component.Sharper cutoffs (larger ∆) prefer softer spectra (larger γ) and a larger rigidity cutoff at the sources, for both hadronic models.Actually, the different parameters in each case combine to produce a similar injection spectral shape at the sources in the most relevant energy range, as it can be seen in Fig. 9 in Appendix B. The deviance obtained is smaller for the milder (∆ = 1) cutoff, although requiring a very hard spectrum, and the deviance grows for increasing ∆, which lead to softer spectral indices.For all the cases considered we obtain that γ H < 1, and in particular for ∆ = 1 one has that γ H < −1.67.For a given value of ∆ the fits using EPOS-LHC lead to smaller deviances than those considering Sibyll 2.3d, while these latter lead to larger values of γ H .The cutoff of the low-energy population often slides to the maximum allowed value in the fit (100 EeV), although the deviance is quite degenerate for values of R L cut larger than 25 EeV since the LE population makes a subdominant contribution at those energies.We report in these cases the range of R L cut leading to a deviance within one unit from the value at the boundary.

Results obtained including the magnetic horizon effect
When including the magnetic horizon effect, two new parameters enter in the fit, the normalized intersource distance X s and the critical rigidity R crit .Note that since for intersource distances much smaller than the diffusion length in the EGMF the spectrum of the particles reaching Earth is not modified (in agreement with the propagation theorem [19]), in the limit X s → 0 the results of the fit should coincide with the ones in the absence of magnetic fields presented above.To gain insight on the dependence of the results with the normalised inter-source distance parameter X s , we first present some results fixing different X s values and fitting the rest of the parameters.We display in Fig. 2 (top panel) the deviances as a function of X s for the different cases considered of non-evolving sources, i.e. for the two hadronic models and for the cutoff functions with ∆ = 1, 2 and 3.One can see that in the case with ∆ = 1 the deviance is degenerate for all X s values, since the fit actually favours the no magnetic field case (i.e.being compatible with a vanishing R crit ), and preferring a negative γ H for all values of X s .For the cases with ∆ = 2 and 3 the fit actually favours the case where the magnetic horizon effect plays a significant role, since smaller deviances are obtained for larger X s , and the deviance becomes almost degenerate for X s ≥ 2. The smallest deviance is obtained for ∆ = 1 (no magnetic horizon), but let us mention that when including the possible impact of systematic experimental uncertainties, the models leading to the smallest deviance may change, as will be discussed in the next Section.In the four lower panels of the figure we show with blue bands the best fitting spectral index of the HE population and with red bands the critical rigidity as a function of X s for the cases with ∆ = 2 and 3, for which the presence of an EGMF is preferred.Here softer associated values of the spectral indices γ H result for increasing X s as a consequence of the magnetic horizon effect.It is interesting to note that larger normalised distances also give rise to smaller best-fitting critical rigidities, since smaller magnetic field strengths are required to produce the spectral suppression at low energies when the sources are farther apart, as it was shown in Fig. 1.We also show in the plots (analytic dashed lines) the approximate dependence of R crit as a function of X s in the region where the magnetic horizon effect plays a relevant role.In the four cases (EPOS-LHC and Sibyll 2.3d with ∆ = 2 or 3), the best fits are obtained for models satisfying Table 2. Parameters of the fit to the flux and composition including the magnetic horizon effect for the EPOS-LHC and Sibyll 2.3d hadronic interaction models and for different cosmological evolutions of the lowenergy component and NE for the high-energy component.Cutoff shapes with ∆ = 1, 2 and 3 are considered.
In Table 2 the best fitting values of the spectral and magnetic suppression parameters, obtained by minimizing also on X s within the range [0, 4], are reported.For EPOS-LHC the minimum deviance corresponds to the ∆ = 1 cutoff case (for which the fit is equivalent to the no magnetic field case), but for Sibyll 2.3d the deviance is smaller for the steeper cutoffs.For ∆ = 2 and 3 a significantly softer spectrum of the HE population with respect to the case with no magnetic field is obtained, arising from the effect of the magnetic horizon which significantly hardens the spectrum of the CRs arriving to the Earth and hence allows for softer spectra at the sources.The spectral indices obtained in this case are larger than unity, and in particular for Sibyll 2.3d with ∆ = 3 a value γ H ≃ 2 is obtained, which is well within the expectations from diffusive shock acceleration.
We also report in Table 2 the results of the fit when considering a scenario involving SFR evolution for the LE component and NE for the HE component (labelled as SFR-NE).This scenario leads to slightly harder LE spectrum, with γ L being smaller by about 0.3, due to the enhanced emission at high redshifts that leads to more significant steepening of the spectrum at the Earth than the NE case due to propagation effects.The rest of the parameters and the deviances are very similar to the NE-NE case.Considering scenarios with SFR evolution for the high-energy population leads always to worse deviances (in agreement with [2]) and thus they will not be further discussed here.
Figures 3 and 4 present the spectrum at Earth and the first two X max moments for EPOS-LHC and Sibyll 2.3d, respectively.The left panels show the results for a cutoff with ∆ = 1, for which a significant magnetic horizon effect is not favoured by the fit.The right plots are for ∆ = 3, where the inclusion of the magnetic suppression leads to a better fit.The case with ∆ = 2 is qualitatively similar to that with ∆ = 3.Although the full X max distribution was fitted, as explained in Section 4.2, we display for illustration of the goodness of the fit the results for the first two X max moments.
Figure 5 shows the differential particle generation rate per logarithmic energy bin of each component at the sources for the different scenarios discussed above.The LE population component, depicted in solid lines, show a soft spectrum for the cases.Note in particular the difference in the spectrum of the HE component (in dotted lines) between the left panels corresponding to ∆ = 1, where there is no magnetic horizon for the best fit and the spectrum at the source is very hard, and the right panels corresponding to ∆ = 3, where the magnetic horizon is significant and the spectrum is much softer.Some features which can be inferred from the results of the fits are: • In all cases the lighter (H and He) nuclei reaching the Earth at energies of few EeV are actually mostly secondaries, arising mainly from the disintegration of N nuclei from the high-energy population.
• There is a significant amount of secondary protons around 1-3 EeV, which leads to a relatively light composition in this energy range.
• The He nuclei are mainly responsible for the instep feature (the change in slope observed at about 15 EeV), whose energy is to some extent related to the energy at which He nuclei get strongly disintegrated by CMB photons but is also affected by the source cutoff suppression.However, in the case of Sibyll 2.3d with ∆ = 3 the N contribution is more significant at the instep and is also partially responsible for this feature.
• Beyond the instep feature, the spectrum is largely dominated by the N component.The suppression above 50 EeV is in part due to the attenuation of the N flux by interactions with the CMB and in part due to the cutoff of the high-energy population, with the cutoff suppression of the N nuclei appearing at an energy ∼ 7/2 = 3.5 times larger than the He one.At the highest energies the Si and Fe elements give the main contribution to the flux.Note that above 10 19.6 eV there is just one integral energy bin of X max , so that the composition information here 10 38 10 38  Flux at Earth (upper panels) and moments of the Xmax distribution (lower panels) for the best fit models, adopting the EPOS-LHC hadronic model.All the scenarios assume NE-NE for the cosmological evolution of the sources.Dotted lines in the spectrum plots represent the flux coming from the primary nuclei while solid lines correspond to the total flux (primary plus secondary) of each mass group.The left panel depict the results for the ∆ = 1 cutoff case, for which the fit prefers the no magnetic horizon solution, while the right panel depict the ∆ = 3 cutoff case, where the best fit has a significant magnetic horizon.10 38 10 38  is quite limited, and improvements in this respect are expected to be obtained with the use of the surface detector information and the recent upgrade of the Observatory [49].
• Regarding the low-energy component, most of the flux at EeV energies arises from the H and N components, with a sizeable contribution also from He in the case of Sibyll 2.3d, and no significant amounts of Si or Fe.
• In the cases where the cutoff rigidity of the low-energy component slides to the limit allowed in the fit (R L cut = 100 EeV) the deviance is actually almost unchanged for values larger than ∼ 25 EeV.These large cutoff values give rise to a subdominant component of H nuclei extending up to very high energies, which could have implications for the production of cosmogenic neutrinos up to EeV energies [2].We also note that often a secondary minimum of the deviance exist with a lower cutoff rigidity R L cut ≃ 3 EeV and a deviance in general larger by a few units, as discussed in [2], although this minimum is actually the global one for Sibyll 2.3d and ∆ = 1.
• The EPOS-LHC hadronic model generally predicts that showers of a given composition and energy are less penetrating (have smaller X max ) than those from the Sibyll 2.3d model.This leads to a lighter inferred composition when the EPOS-LHC model is considered, as is apparent in the fact that for this model there is an enhanced contribution of H observed below the ankle, the He contribution is more prominent around the instep and a slightly smaller Si contribution is present in the suppression region.
We report in Table 3 the results of the minimization for the full set of fitted parameters, including their statistical uncertainties, for the scenarios shown in Figs. 3 and 4. They correspond to one case where the magnetic horizon effect plays an important role (∆ = 3) and another where it does not (∆ = 1), both for the EPOS-LHC and Sibyll 2.3d hadronic models.We also quote the present day emissivities above a threshold energy E th = 10 17.8 eV, L a 44 ≡ A L a A /(10 44 erg Mpc −3 yr −1 ), where It is worth noting that the emissivities inferred for the HE population in the scenarios where the magnetic horizon effect plays a significant role are larger than in the absence of EGMF, since in the first case there is a sizeable emission at low energies which doesn't manage to arrive to the Earth.

Impact of systematic uncertainties on the energy and X max calibrations
We explore in this section the effect that the experimental systematic uncertainties on the energy scale and X max calibration have on the fit results.For the energy scale, an energy independent uncertainty ∆E/E = ±14% is considered in the whole energy range analysed [48].The systematic uncertainties on the measured X max values are asymmetric and slightly energy-dependent, ranging from 6 to 9 g cm −2 [8].We quantify the effects of these uncertainties by shifting the measured energies and the inferred X max values by one systematic standard deviation up and down, and performing the fit again for the nine possible combinations of shifted and unshifted data sets.Figure 6 displays the resulting total deviance and the HE population spectral index obtained for the different cutoffs and hadronic models considered.In general the deviance improves for a positive shift in energy and/or a negative shift in X max .On the other hand, a positive shift in X max leads to a deviance larger by more than 100 units, specially worsening the fit to the composition data.For completeness we include in Appendix C the results for the fits performed under similar shifts of energy scale and X max in the absence of EGMF, where similar trends are observed, although generally the deviances are larger.Since Sibyll 2.3d with ∆ = 3 leads to the smaller deviance when performing a positive 1 σ shift in the energy and a negative 1 σ shift in X max , with a γ H ≃ 2, we display in Table 3. Parameters of the fit to the flux and composition for the scenarios shown in Figs. 3 and 4. We include for the EPOS-LHC and Sibyll 2.3d hadronic interaction models a scenario with ∆ = 1 (for which the magnetic horizon effect is not relevant), as well as a scenario with ∆ = 3 and with EGMF for NR of source luminosities.Quoted are the fitted Xs and Rcrit, the spectrum shape parameters and element fractions for the two populations as well as the source emissivities above 10 the rigidity cutoff of the high-energy component moves to larger values for a positive energy shift (and to smaller values for a negative shift), while γ H ≃ 2 holds in all cases.The low energy spectral parameters remain practically unchanged for most of the shifts.On the other side, a positive shift in  impact on the combined fit results when taking into account the magnetic horizon effect.We included this effect only for the high-energy component, given its assumed lower source density, considering that the associated effects in the denser low-energy component happen at energies below those considered in the analysis.In particular, if the intersource separation of the low-energy component were an order of magnitude smaller than that of the high-energy one, the magnetic suppression effects on it would be very mild above 1 EeV.Let us note that taking into account that the inferred emissivity per unit volume and time for the low and high energy populations are of comparable magnitude, one expects in that case that the typical individual luminosities of the much more abundant LE sources be fainter by at least three orders of magnitude with respect to those of the HE sources.We found that, in order for the magnetic horizon effect to play a relevant role, the normalised intersource distance should satisfy X s > 1, and hence L coh 100 kpc . (5.1) Larger values of X s require smaller critical rigidities, and the deviance is almost degenerate for X s > 2.
For the scenarios where the magnetic horizon effect is responsible for the hardness of the inferred HE spectra we obtained that, allowing for the possible systematic shifts on the measurements that were discussed in the previous section, the approximate relation X s R crit ≃ 5 to 10 EeV should hold, where these quantities are related through depends on the four parameters a, b, α and β, which in turn are functions of the spectral index γ, the source evolution and the kind of nuclei considered.The value of the parameters, obtained in [45], are reported in Table 5.
Primaries refer to nuclei that reach Earth as part of the same mass group as that in which they were injected.Secondary protons are H nuclei that were produced via photodisintegration processes.Intermediate secondary nuclei refers to secondary nuclei, produced in photodisintegrations, belonging to a lighter mass group than the primary nuclei.

1
Results obtained in the absence of magnetic fields 7 4.4.2Results obtained including the magnetic horizon effect 7 4.5 Impact of systematic uncertainties on the energy and X max calibrations 12

Figure 1 .
Figure 1.Magnetic suppression factor as a function of energy over critical energy for three values of the normalized intersource distance Xs.The left panel assumes a no evolution scenario while the right panel assumes a star formation rate one.Results are shown for primaries (solid lines), secondary protons (dashed lines) and intermediate secondary nuclei (dotted lines).

Figure 2 .
Figure 2. The top panel displays the deviances as a function of Xs for the different scenarios of the NE-NE case (i.e.NE for the LE and HE components).Lower panels display the values of Rcrit and γH as a function of Xs for the two hadronic models and for ∆ = 2 (left panels) or ∆ = 3 (right panels).Also shown are analytic fits to Rcrit as a function of Xs.

Figure 3 .
Figure 3. Flux at Earth (upper panels) and moments of the Xmax distribution (lower panels) for the best fit models, adopting the EPOS-LHC hadronic model.All the scenarios assume NE-NE for the cosmological evolution of the sources.Dotted lines in the spectrum plots represent the flux coming from the primary nuclei while solid lines correspond to the total flux (primary plus secondary) of each mass group.The left panel depict the results for the ∆ = 1 cutoff case, for which the fit prefers the no magnetic horizon solution, while the right panel depict the ∆ = 3 cutoff case, where the best fit has a significant magnetic horizon.

. 5 )Figure 5 .
Figure 5. Differential injection rate of particles at the source per logarithmic energy bin for H (red), He (gray), N (green), Si (cyan) and Fe (blue) for the scenarios considered in Figs.3 and 4: the EPOS-LHC and Sibyll 2.3d hadronic interaction models and the ∆ = 1 and ∆ = 3 cutoff shapes for NE of the source luminosities.

Table 4 a
more in-depth exploration of the effect of the systematic shifts on the fit parameters for this scenario.As expected,

Table 4 .
17.8eV expressed in units of 10 44 erg Mpc −3 yr −1 .The bottom row indicates the associated deviances and the number N of data points considered.Effects of systematic uncertainties: results for shifts of ±σsys in the energy and Xmax scales.

Table 5 .
Parameters of the fit for the magnetic suppression factor G, using eq.(3.2), for primary nuclei, secondary protons and intermediate mass secondary nuclei, for both the NE and SFR scenarios.