Observational constraints on cosmic-ray escape from ultra-high energy accelerators

The sources of ultra-high energy cosmic rays (UHECRs) remain unknown. Their identification relies primarily on capturing in the UHECR arrival directions a pattern suggestive in an evident way of a class of astrophysical objects. Such a capture is still eluding our grasp, but some recent observations have allowed broad statements to be drawn. An anisotropy at large scales has been revealed above ${\simeq}\>8\>$EeV (Aab et al., 2017a), the contrast and the direction of which are consistent with expectations drawn from sources distributed in a similar manner to the extragalactic matter (Aab et al., 2017a, 2018a). At higher energies (${\gtrsim}\>40\>$EeV), a correlation between UHECR arrival directions and the flux patterns of massive, star-forming or active galaxies within 200 Mpc provide evidence for anisotropy (Aab et al., 2018b; Biteau & Pierre Auger Collaboration, 2021). Overall, these observations suggest that UHECRs are predominantly of extragalactic origin at least above the so-called ankle energy at ${\simeq}\>5\>$EeV.

The energy spectrum and chemical composition of UHECRs observed on Earth result from the emission processes at play, which encompass acceleration mechanisms, losses and escape from the source environments as well as propagation effects. Differently from, and complementary to anisotropies, these two observables provide constraints helping infer the properties of the acceleration processes, the energetics of the sources, and the abundances of elements in the source environments. Following this principle, several studies have laid the foundations for a generic scenario that broadly reproduces the observations (Allard et al., 2008; Aloisio et al., 2014; Taylor et al., 2015; Aab et al., 2017b; Zhang et al., 2018; Guido & Pierre Auger Collaboration, 2021). The intensity of the individual nuclear components at the sources, generally considered as stationary and uniformly distributed in a comoving volume, is assumed to drop off at the same magnetic rigidity so as to explain the gradual increase with energy of mass number $A$ observed on Earth (Abraham et al., 2010; Aab et al., 2014a; Watson, 2022). This is consistent with the basic expectation that electromagnetic processes accelerate particles up to a maximum energy proportional to their electric charge $Z$. Most notably, the abundance of nuclear elements is found to be dominated by intermediate-mass ones, ranging from He to Si, accelerated to $E^{Z}_{\mathrm{max}}\simeq 5Z\,\text{EeV}$ and escaping from the source environments with a very hard spectral index $\gamma$ (that characterizes the emission spectrum at the sources as $E^{-\gamma}$), which, depending on the systematic uncertainties affecting some of the necessary modeling, ranges between $\gamma\simeq-1$ and $\gamma\simeq 1$ (Aab et al., 2017b).

Reaching ultra-high energies can be most easily achieved by first-order Fermi shock acceleration, where particles are energized by bouncing back and forth across moving magnetic fields (see e.g. Sironi et al., 2015, for a review). While the typical spectral index for this mechanism ($\gamma\gtrsim 2$) appears to deviate significantly from the hard values favored from the data, the interplay between the escape from the sources and the acceleration mechanisms must be taken into account to correctly apprehend the emission spectrum (e.g. Fujita et al., 2009), in particular when escape is influenced by in-source interactions (Unger et al., 2015; Globus et al., 2015; Biehl et al., 2018; Zhang et al., 2018; Fang & Murase, 2018; Supanitsky et al., 2018; Boncioli et al., 2019; Muzio et al., 2022). The cosmic-ray luminosity of the source candidates is governed, among other things, by the levels of radiation and magnetic-field densities. While gaining energy, UHECRs can interact with this radiation or escape from the magnetized zone. Thus, the energy spectrum of the ejected particles as well as the amount of ejected nuclei may differ strongly from those injected into the electromagnetic field. This has two important consequences: a) the ejected spectrum of the charged nuclei can be much harder than that injected, due to the escape mechanism and in particular the behavior of the nuclei-photon cross section at high energies; b) the interactions can produce a copious flux of secondary neutrons of energy $E_{n}=E/A$. These neutrons can escape freely from the magnetic confinement zones, with an ejection spectrum much softer than that of nuclei, to decay into protons on their way to the Earth after an average travel distance of $9.1\>\text{kpc}\times(E_{n}/1\>\text{EeV})$.

In this paper, we derive constraints on the characteristics of the sources by searching for a robust statistical agreement between data and predictions based on realistic astrophysical ingredients and on a phenomelogical model incorporating in-source interactions in an effective way. The agreement is sought for both spectral and mass composition data, while only a qualitative overview of the latter has been performed so far in studies incorporating in-source interactions. To do this, we consider the data obtained at the Pierre Auger Observatory beyond 0.63 EeV ($10^{17.8}~{}$eV). However, and differently from the approaches adopted elsewhere, only the proton spectrum is used in the energy range between 0.63 EeV and 5 EeV. This approach allows us to reconstruct the extragalactic component without resorting to any introduction of ad hoc nuclear components to model what is generally assumed to be the upper end of the galactic component so as to reproduce the average mass composition as a function of energy. Such a modeling, suffering from both observational systematics and theoretical unknowns, blurs the reconstruction of the extragalactic proton component with biases inherent to the choices made. In the same spirit, we do not attempt to sew together the Auger observations with those made in an energy range half-a-decade lower with the KASCADE-Grande detector mixing indistinctly protons and Helium (Apel et al., 2012), always with the aim of reconstructing the proton component as precisely as possible.

With this approach, we show that a scenario incorporating in-source interactions can reproduce the data and that accounting for the local overdensity of galaxies within tens of Mpc further improves the goodness of fit with respect to the generally assumed uniform distribution of sources in a comoving volume. Quantitatively, we show that the analysis of data from the Pierre Auger Observatory significantly favors an effective spectral index of protons ejected from the sources that is softer than generally obtained or assumed, somehow pushing the in-source interaction mechanism to the extreme. At the same time, we point out that this spectral index approaches more “canonical” values provided that all protons are secondaries of the interactions, i.e. that the environment of the sources is devoid of protons to be accelerated.

This paper is organized as follows. In Section 2, we motivate the model used to probe those effects and set up a benchmark scenario among the ingredients needed to propagate the particles from their sources to the Earth. The data and fitting method are presented in Section 3. Results obtained with the benchmark scenario are given in Section 4, where we explore as well the impact of changing the ingredients that are affected by systematic uncertainties. Finally, the significance of the results is discussed in Section 5.

The benchmark astrophysical model used here is inspired by that of Aab et al. (2017b) for its main features. The non-thermal processes responsible for accelerating the different particles are modeled through power-law spectra as long as energies are sufficiently below the $Z$-dependent maximal acceleration energy $E^{Z}_{\mathrm{max}}=ZE_{\mathrm{max}}$, while, in the absence of any firmly-established prescription from theory, an exponential suppression is used for modeling the upper end of the acceleration process. The sources are assumed to accelerate different amounts of nuclei represented by five stable ones: hydrogen (¹H), helium (⁴He), nitrogen (¹⁴N), silicon (²⁸Si) and iron (⁵⁶Fe). The unavoidable fluctuations of the characteristics of each individual source are neglected, considering all sources as identical. This is a simplification that renders the number of free parameters manageable in a fitting procedure. The fitted parameters should be interpreted as effective ones, describing the average spectrum of the source population.

Modeling the ejection spectra in terms of power laws suppressed at the highest energies is an approximation aimed at keeping, once again, the number of free parameters to a minimum. However, detailed studies of the nuclear cascade developing in sources from photodissociation with different levels of radiation densities, such as that presented in Biehl et al. (2018), show that such a simplification holds in the energy range of interest for this study. This is also true for the softer flux of secondaries neutrons produced during the development of the cascade, which dominates over that of ejected protons. We generically model the ejection rate per comoving volume unit and per energy unit of nucleons as

with $Z_{\rm p}=1$, that is a single reference ejection rate $q_{0{\rm p}}$, spectral index $\gamma_{\rm p}$ and suppression function $f_{\mathrm{supp}}(E,Z_{\rm p})$ for both escaping protons and protons from neutron decay. Here, $E_{0}$ is arbitrarily set to 1 EeV. The ejection rate of nuclei with mass number ${A}_{i}$ is also generically modeled as

that is a single spectral index $\gamma_{A}$ and four independent reference ejection rates $q_{0{A}_{i}}$ for helium, nitrogen, silicon and iron. As in the reference scenario of Aab et al. (2017b), the suppression function adopted both for nucleons and nuclei is taken as

The maximum acceleration energy is assumed to be proportional to the electric charge of each element, $E^{Z}_{\mathrm{max}}=ZE_{\mathrm{max}}$, with a single free parameter $E_{\mathrm{max}}$ shared by the five species.

In this approach, the ejection rate for protons accounts for both the accelerated ones up to $E_{\mathrm{max}}$ and those produced by the escaping neutrons, inheriting hence an energy of the order of $E/A$ from the nuclei of energy $E$. The maximum energy that this population of secondary protons can reach is then $E^{Z}_{\mathrm{max}}/A\simeq E_{\mathrm{max}}/2$, which is not reflected in the rigidity-acceleration scheme modeled by Equation (3). The extreme case in which all ejected protons would be photodissociation by-products can be tested by replacing $E_{\mathrm{max}}$ by $E_{\mathrm{max}}/2$ in $q_{\mathrm{p}}(E)$. We will explore this extreme case in Section 4. Further characterization of the balance between the population of accelerated protons from the initial abundance in the source environment and that of secondaries from nuclear cascades calls for modeling of specific sources and environments, which is beyond the scope of this paper.

The differential energy production rate per comoving volume unit of the sources, which is directly related to their differential luminosity, is then $\ell_{j}(E,z)=E^{2}q_{j}(E)S(z)$ for each species $j$, where $S(z)$ describes the redshift evolution of the UHECR luminosity density. The quantity $\ell_{j}(E,z)$ is hereafter called differential energy production rate for convenience. On the other hand, the bolometric energy production rate per comoving volume unit at redshift $z$ is obtained as $\mathcal{L}_{j}(E,z)=S(z)\int_{E}^{\infty}\mathrm{d}E^{\prime}E^{\prime}q_{j}(E^{\prime})$. We hereby report its average value in a volume spanning $z_{\text{min}}-z_{\text{max}}$ as $\bar{\mathcal{L}}_{j}(E)=\int_{z_{\text{min}}}^{z_{\text{max}}}\mathrm{d}z\left|\frac{\mathrm{d}t}{\mathrm{d}z}\right|\mathcal{L}_{j}(E,z)/\int_{z_{\text{min}}}^{z_{\text{max}}}\mathrm{d}z\left|\frac{\mathrm{d}t}{\mathrm{d}z}\right|$, where $t(z)$ is the lookback time.

The evolution of the UHECR luminosity density is taken as being traced by the density of baryonic matter over cosmic time. The latter is assumed to follow the density of stellar mass, which is fairly approximated by a constant out to redshift $z=1$ (e.g. Madau & Dickinson, 2014). As in the benchmark scenario of Aab et al. (2017b), such a constant evolution is assumed to hold to first approximation out to $z_{\text{max}}=2.5$. The local Universe presents however an overdensity because, like most galaxies, the Milky Way belongs to a group of galaxies, itself embedded in the Local Sheet (McCall, 2014). We adopt here the overdensity correction factor inferred by Condon et al. (2019) expressed as a function of the distance $r$ and effective up to 30 Mpc,

with $r_{0}=5.4\>$Mpc and $\alpha=1.66$. The evolution of the UHECR luminosity density is then described as $S(z)=\delta\rho/\bar{\rho}$, with $z_{\text{max}}=2.5$ and $z_{\text{min}}=2\times 10^{-4}$. The minimum redshift, which corresponds to the limit of the Local Group of galaxies ($r\simeq 1\>$Mpc), prevents any divergence in Equation (4) and effectively removes very-closeby galaxies that would otherwise dominate the UHECR sky, at odds with observations.

Extragalactic magnetic fields are poorly known, except for upper limits at the nG level from rotation measures (Pshirkov et al., 2016; O’Sullivan et al., 2020) or down to tens of pG for magnetic fields of primordial origin that would affect CMB anisotropies (Jedamzik & Saveliev, 2019). Lower limits at the fG level have also been derived from the non-observation in the GeV range of gamma-ray cascades from TeV blazars (Neronov & Vovk, 2010; Tavecchio et al., 2010; Ackermann et al., 2018). These non-observations have been interpreted as the deflection of the cascade flux into a broadened beam weakening the point-like image, noting though that the lower limits could be alleviated if plasma instabilities dominated the cooling rate of particles in the cascade (see Alves Batista & Saveliev, 2021; Biteau & Meyer, 2022, for recent reviews). While field values saturating the limits from rotation measures might impact the results by reducing the UHECR horizon (Mollerach & Roulet, 2013; Wittkowski, 2018), we assume here that magnetic fields in cosmic voids are at the level of 0.1 nG or less, as suggested by cosmological magneto-hydrodynamical simulations (Vazza et al., 2017). In this case, they do not have sizeable effects on the propagation of the particles, which can therefore be considered one-dimensional.

The all-particle energy spectrum $J(E)$ observed at present time thus results from the integration of the contribution of all sources over lookback time, the role of which is played by redshift:

Here, the relationship between cosmic time and redshift follows from the concordance model used in cosmology, $(\mathrm{d}t/\mathrm{d}z)^{-1}=-H_{0}(1+z)\sqrt{\Omega_{\mathrm{m}}(1+z)^{3}+\Omega_{\Lambda}}$, where $H_{0}=70\>$km s^-1 Mpc^-1 is the Hubble constant at present time, $\Omega_{\mathrm{m}}\simeq 0.3$ is the density of matter (baryonic and dark matter) and $\Omega_{\Lambda}\simeq 0.7$ is the dark-energy density. The energy losses and spallation processes are described by $\eta_{AA^{\prime}}(E,E^{\prime},z)$, which is the fraction of particles detected on Earth with energy $E$ and mass number $A$ from parent particles emitted by the sources with energies $E^{\prime}>E$ and mass numbers $A^{\prime}>A$. In practice, for a given source with redshift $z_{0}$ emitting a nuclear species $A_{0}$ at energy $E_{0}$, the corresponding $\eta_{A_{0}A^{\prime}}(E_{0},E^{\prime},z_{0})$ function is tabulated in bins of $A^{\prime}$ and $E^{\prime}$ by propagating a large number of emitted particles (O($10^{7}$) particles) using the SimProp package (Aloisio et al., 2012, 2015, 2016). By repeating the simulations for different values of $z_{0}$, $A_{0}$ and $E_{0}$, the whole $\eta_{AA^{\prime}}(E,E^{\prime},z)$ function is tabulated as a 5D histogram, providing hence the fractions sought for. The relevant processes accounted for in SimProp are pair production, photo-pion production and photodissociation off the photon fields of interest, which are those from the cosmic-microwave-background (CMB) radiation, and the infrared photons from the extragalactic background light (EBL) that comprises the radiation produced in the Universe since the formation of the first stars. The CMB radiation is characterized as a black-body spectrum with a redshift-dependent temperature $T(z)=T_{0}(1+z)$, with $T_{0}=2.725\>$K. The EBL is less precisely known, especially in the far infrared and at high redshifts. For the benchmark scenario explored below, we use the model of Gilmore et al. (2012) that is consistent with the minimal intensity level from galaxy counts and with indirect measurements from absorption of gamma-rays at multi-TeV energies (e.g. Pueschel & Biteau, 2021, for a recent review). Alternative EBL models are studied as sources of systematic uncertainties. Pair production is a very well-known process that results in a small fractional energy loss above a threshold of ${\simeq}\>0.5A\epsilon_{\mathrm{meV}}\>$EeV, with $\epsilon_{\mathrm{meV}}$ the energy of the photon targets in meV. Photo-pion cross sections have been measured in accelerator-based experiments and are well reproduced by various event-generator codes. The corresponding fractional energy loss is quite important above a threshold of ${\simeq}\>30A\>$EeV, causing the GZK effect. Photodissociation cross sections for nuclei, on the other hand, are less known, especially for exclusive channels in which charged fragments are ejected. The fractional energy losses used here, which also shape the suppression of the spectrum, rely on phenomenological approaches, using the TALYS model (Koning et al., 2005; Koning & Rochman, 2012) for the benchmark scenario. We also evaluate alternative models as sources of systematic uncertainties. In any event, because the binding energies are of order of a few MeV for all nuclei of interest, the center-of-mass energy depends only on the Lorentz factor of the nuclei. For heavy nuclei like iron, the photodissociation thus occurs on higher-energy photons – infrared background photons with energies about an order of magnitude higher than the typical CMB photon energy. Since these photons are less abundant than the CMB ones, the attenuation of heavy nuclei is relatively slow and is comparable to that of protons. On the other hand, the Lorentz factor is high enough for light nuclei to photodissociate on the CMB background and their attenuation is much faster than that of protons. Finally, the adiabatic losses due to the expansion of the Universe are included in the energy-loss rate of the particles as $-(1/E)(\mathrm{d}E/\mathrm{d}t)=-((1+z)\mathrm{d}t/\mathrm{d}z)^{-1}$.

With these ingredients, Equation (5) can be used to evaluate the all-particle flux from the various contributions of each individual nuclear component on the condition to assign values to the five ejection rates $q_{0A_{i}}$, the two indices $\gamma_{\rm p}$ and $\gamma_{A}$, and the maximum energy $E_{\rm max}$. These eight parameters are fitted to the data in the way explained in the next section.

The expected spectrum modeled by Equation (5) depends on several unknown parameters characterizing properties of the acceleration processes, of the source environments and of the source energetics. The observed energy spectrum and mass composition can, as discussed in Section 1, provide constraints helping in inferring these parameters.

The data we use hereafter are those obtained at the Pierre Auger Observatory, located in the province of Mendoza (Argentina) and providing the largest exposure to date to UHECRs (Aab et al., 2015). The Observatory is a hybrid system that detects the extensive air showers induced in the atmosphere subsequent to the collisions of UHECRs with nitrogen and oxygen molecules. We shall make use on the one hand of the all-particle energy spectrum inferred from these data (Abreu et al., 2021), and on the other hand of the distributions of the slant depth of maximum of shower development ($X_{\mathrm{max}}$), which is a proxy, the best up to date, of the primary mass of the particles (Aab et al., 2014a; Bellido & Pierre Auger Collaboration, 2018). With caution over the hadronic-interaction generators used to model the development of the showers, the $X_{\mathrm{max}}$ distributions allow for inference of the energy-dependent mass composition on a statistical basis. Two hadronic-interaction generators are considered here, namely EPOS-LHC (Pierog et al., 2015) for the benchmark scenario and Sibyll2.3c (Riehn et al., 2017) as an alternative, which are those up-to-date that best describe the data. Hence, the various mass components $J_{j}(E)$ can be derived by combining the all-particle energy spectrum and the abundances of the different elements as a function of energy. The energy threshold considered here is the nominal one of the detection mode of $X_{\mathrm{max}}$ at ultra-high energies, namely ${\simeq}\>0.63\>$EeV ($10^{17.8}\>$eV). Incidentally, this threshold allows us to explore the energy range of the ankle feature.

From a set of proposed parameters $\mathbf{\Theta}$, the best match between the observed spectrum and that expected from Equation (5) is obtained through a Gaussian likelihood fit. Under the assumption that the transition to extragalactic UHECRs is already completed above 5 EeV,¹¹1We checked that a moderate increase in the value at which the transition to extragalactic UHECRs is assumed to be completed has no significant impact on the results. the corresponding likelihood term, $L_{J}$, is calculated above that threshold. The model is fitted to the $X_{\mathrm{max}}$ data, following Aab et al. (2017b), using the multinomial distribution that describes how likely it is to observe, in each energy bin $m$, $k_{mx}$ events out of $n_{m}$ with probability $p_{mx}$ in each $X_{\mathrm{max}}$ bin. The probabilities $p_{mx}$ are obtained by using generators of hadronic interactions to model the $X_{\mathrm{max}}$ distributions expressed in terms of the proposed parameters $\mathbf{\Theta}$. The corresponding likelihood term, $L_{X_{\mathrm{max}}}$, is built as the product over energy bins above 5 EeV. The last contribution to the likelihood stems from the sub-ankle proton component, $L_{J_{\mathrm{p}}}$, obtained by weighting the all-particle spectrum below 5 EeV with the proton abundance $f_{\mathrm{p}}(E)$ from the $X_{\mathrm{max}}$ distributions. The statistical uncertainties in $\hat{J}_{\mathrm{p}}(E)$ are dominated by those in $f_{\mathrm{p}}(E)$ and are accounted for through a Gaussian likelihood fit.

The model likelihood is therefore given by $L=L_{J}L_{X_{\mathrm{max}}}L_{J_{\mathrm{p}}}$. The goodness-of-fit is assessed with a deviance, $D$, defined as the negative log-likelihood ratio of a given model and the saturated model that perfectly describes the data:

The three different contributions are referred to as $D_{J}$, $D_{X_{\mathrm{max}}}$ and $D_{J_{\mathrm{p}}}$ respectively.

Fitting the model to the data within the framework of the benchmark scenario described in Section 2, with free spectral indices for both nucleons and nuclei ($\gamma_{\rm p}\neq\gamma_{A}$), leads to the parameters and deviance given in Table 1 using EPOS-LHC to interpret $X_{\mathrm{max}}$ data (results obtained using Sibyll2.3c are given in Appendix C for completeness). For reference, Table 1 also provides the results obtained under the assumption of a proton component dominated by neutron escape (proton maximum energy of $E_{\rm max}/2$), of no local overdensity (widely-used uniform distribution) and of a shared spectral index across the five species ($\gamma_{\rm p}=\gamma_{A}$).

Similarly to that found in Aab et al. (2017a), the value of $E_{\mathrm{max}}$, determined by the drop in the nuclear components at $E^{Z}_{\mathrm{max}}$, implies that the suppression of the spectrum is due to the combination of the cut-off energy at the sources for the heavier nuclei and the energy losses en route. The spectral index of the nuclei, $\gamma_{A}$, is in turn determined by the increase of the average mass with energy, which is almost monoelemental, so as to reproduce the $X_{\mathrm{max}}$ distributions as well as possible. The solution provided by the fit therefore consists in imposing a hard index for nuclei so that the contribution of each element mixes as little as possible: high-energy suppression imposed by the cut-off beyond $E^{Z}_{\mathrm{max}}$ and low-energy suppression via the hard index $\gamma_{A}$. However, this phenomenon does not apply to protons, which are present in an energy range where a mixture of elements is required. The best-fit value of $\gamma_{\mathrm{p}}$ is much softer than that of $\gamma_{\mathrm{A}}$. Note that the introduction of $\gamma_{\mathrm{p}}$ significantly improves the fit of the data down to 0.63 EeV, with a total deviance $D=236.8$ compared to $862.7$ in the case of $\gamma_{\mathrm{p}}=\gamma_{\mathrm{A}}$ so that the introduction of this extra free parameter is amply justified. On the other hand, the introduction of the local overdensity to trace the source distribution also provides a substantial improvement of the deviance, although with moderate impact on the best-fit parameters.

The balance regulating the intensity of each component is reported in terms of the energy production rates $\bar{\mathcal{L}}$ above 0.63 EeV. The solution is illustrated in Figure 1, where the contributions of each nuclear component to the observed energy flux and energy density are displayed. In Figure 2, the energy production rates required at the sources to fuel the observed energy flux are shown as a function of energy for the different primary mass groups. While the contribution of nuclei peaks at most to ${\simeq}\>5\times 10^{44}\>{\rm erg}\>{\rm Mpc}^{-3}\>{\rm yr}^{-1}\>{\rm dex}^{-1}$, that of protons is increasing up to ${\simeq}\>10^{45}\>{\rm erg}\>{\rm Mpc}^{-3}\>{\rm yr}^{-1}\>{\rm dex}^{-1}$ when going down in energy. Extrapolations of the results below 0.63 EeV are however hazardous, as the functional shape used in Equation (1) may not hold anymore depending on the specifics of the source environments that govern the nuclear cascade. Integrated above 0.63 EeV, the total energy production rate is found to be $(10.8\pm 0.4)\times 10^{44}\>{\rm erg}\>{\rm Mpc}^{-3}\>{\rm yr}^{-1}$. For completeness, the $X_{\rm max}$ distributions used in this work together with the best-fit models obtained within the benchmark scenario explored here are shown in Appendix A (Figure 3).

The reduced deviance is decomposed into three terms, according to Equation (6). As detailed in Appendix B, the spectral sector leads to an acceptable fit ($D_{J}+D_{J_{\mathrm{p}}}=37.3$ for $N_{J}+N_{J_{\mathrm{p}}}=24$ points). As in Aab et al. (2017a), the $X_{\mathrm{max}}$ sector is more difficult to fit as evidenced by the value $D_{X_{\mathrm{max}}}=199.5$ for $N_{X_{\mathrm{max}}}=109$ points. The large value of the deviance on $X_{\mathrm{max}}$ reflects the difficulty for hadronic-interaction generators to reproduce the data, as illustrated in Appendix A. In general, the benchmark scenario fits the data with similar performances to those obtained by considering only data beyond 5 EeV. Considering the variation of Equation (1) that consists in substituting $E_{\mathrm{max}}$ for $E_{\mathrm{max}}/2$ yields the results reported in the fourth line of Table 2. It is interesting that the main changes concern the spectral index of protons, which becomes $\gamma_{\mathrm{p}}\simeq 2.5$. The astrophysical consequences of these results will be discussed in Section 5, after the sources of systematic uncertainties have been discussed below.

The inferred mass-dependent energy production rates, tracing directly the abundances of elements escaping from the source environments, rely primarily on the $X_{\mathrm{max}}$ distributions. The systematic uncertainties in $X_{\mathrm{max}}$, which are slightly energy-dependent and range from 6 to 9 g cm^-2 (Aab et al., 2014b), are thus expected to impact the results. A reduction of the $X_{\mathrm{max}}$ scale by $1\>\sigma$, which results in an overall heavier composition on Earth, has as main impact a decrease of the proton and helium contributions by respectively ${\simeq}\>15\%$ and ${\simeq}\>35\%$, with an increase of the spectral index $\gamma_{A}$ of $+0.7$. This configuration leads to a deviance improvement of ${\simeq}\>16$ units. On the other hand, an increase of the $X_{\mathrm{max}}$ scale by $1\>\sigma$ deteriorates the quality of the fit in an untenable way (more than 100 units of deviance), as the overall lighter composition implied is then in tension with the width of the $X_{\mathrm{max}}$ distributions. By contrast, the inferred parameters and quality of the fits are mildly altered by the systematic uncertainties affecting the energy spectrum, which are dominated by those in the energy scale ($\Delta E/E=14\%$). This is because a change in $E_{\mathrm{max}}^{Z}$ can be reproduced by slightly different balances regulating the intensity of each nuclear component due to the dependence of the photodissociation threshold with the mass number $A$ of the nuclei.

The interpretation of the $X_{\mathrm{max}}$ distributions heavily relies on the hadronic-interaction generator used in the simulations of extensive air showers. Based on the generator Sibyll2.3, increased (decreased) fractions of intermediate nuclei (protons) are inferred at all energies (Bellido & Pierre Auger Collaboration, 2018). An overall heavier composition is thus implied, in a manner similar to that obtained with the analysis based on EPOS-LHC after a $-1\>\sigma$ shift of the $X_{\mathrm{max}}$ scale is applied. The results reported in Appendix C are thus similar to those discussed above when applying the $-1\>\sigma$ shift of the $X_{\mathrm{max}}$ scale.

The range of $E_{\mathrm{max}}^{Z}$ makes the role of the EBL more important than that of the CMB to control the energy losses by photodissociation of the nuclei. Consequently, the uncertainties affecting both the far-infrared intensity of the EBL and the partial cross sections of channels in which $\alpha$ particles are ejected can alter the results. Considering the EBL model of Dominguez et al. (2011) instead of Gilmore et al. (2012) results in an increased far-infrared density of photons, which enhances the intensity of secondary protons from photodissociations en route. Most notably, this enhancement is compensated in the fit to the data through an energy production rate of He (CNO) elements that is increased by ${\simeq}\>25\%$ (${\simeq}\>45\%$) while that of protons is decreased by ${\simeq}\>30\%$. Also, the spectral index of protons is increased by $+1$ for a global deviance similar to that of the benchmark scenario. On the other hand, cross sections for photodissociation from Puget et al. (1976); Stecker & Salamon (1999), which neglect $\alpha$-particle production, impact mainly the balance of the energy production rate of He (CNO) [Si], changed by ${\simeq}\>+70\%~{}(-45\%)~{}[+20\%]$, at the cost of a degraded deviance by 12 units.

Finally, we studied the impact of the redshift evolution of the UHECR luminosity density, considered as flat on large-spatial scales in the benchmark scenario. Alternatively to a flat evolution, we tested a scenario where the evolution strictly follows the stellar-mass distribution on large spatial scales (Madau & Dickinson, 2014). The only notable difference with the benchmark scenario is a decrease by 20% of the emissivity of all components, which can be understood from the lower average source distance resulting from the decrease of stellar-mass density with redshift at $z>1$. Another attractive scenario is to assume an evolution proportional to the redshift-dependent star-formation rate. Using the corresponding $S(z)$ function derived in Madau & Dickinson (2014), the most notable differences are an increase of the overall energy production rate by a factor $2.5$ and a change in the balance of H (He) [CNO] {Si} by ${\simeq}\>+40\%(+90\%)[+350\%]\{+25\%\}$. In such a scenario, the emissivity would be dominated by a soft proton component and a hard CNO component, with lower mass nuclei partly generated by photodissociation en route from the cosmic-noon epoch at $z>1$. The quality of the fit is almost unchanged for a source evolution following that of stellar mass or star formation compared to that of the benchmark scenario.

The results presented in the previous section direct the interpretation of the origin of protons below the ankle energy and of the hard spectra of nuclei above this energy. In this picture, the component of protons is of extragalactic origin well below the ankle energy and exponentially suppressed above it, while heavier nuclei steadily take over to the highest energies through a rigidity-dependent maximum-energy scenario. That the protons do not get suppressed when going down in energy at the same fast rate as the heavier nuclei is modeled by a softer spectral index for this population of low-charge primaries. Interestingly, such a behavior qualitatively fits with scenarios of in-source interactions in which copious fluxes of neutrons, which are produced while accelerated charged particles interact with the bath of photons permeating the sources, escape freely the electromagnetic fields.

This work was made possible by with the support of the Institut Pascal at Université Paris-Saclay during the Paris-Saclay Astroparticle Symposium 2021, with the support of the P2IO Laboratory of Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of Université Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master projet UCMN and EuCAPT ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038). SM, JB, AC & OD gratefully acknowledge funding from ANR via the grant MultI-messenger probe of Cosmic Ray Origins (MICRO), ANR-20-CE92-0052. The authors would like to thank colleagues from the Pierre Auger Collaboration for fruitful suggestions, in particular J. M. Gonzalez, S. Mollerach and E. Roulet for pointing out differences in the normalization of the overdensity parametrization in its cumulated and differential forms, S. Petrera for providing the latest parameterizations of $X_{\rm max}$ distributions using various hadronic interaction models, and R. Clay and A. Di Matteo for their careful reading and suggestions on the text.