SimProp v2r4: Monte Carlo simulation code for UHECR propagation

TL;DR

SimProp v2r4 is an improved Monte Carlo simulation tool for ultra-high energy cosmic ray propagation, now including secondary particle production, enabling better interpretation of observational data.

Contribution

The paper introduces SimProp v2r4 with enhanced capabilities and reduced computation time, incorporating secondary cosmogenic particles for more comprehensive cosmic ray modeling.

Findings

Reduced computation time compared to previous versions

Includes secondary particle production during cosmic ray propagation

Provides insights into cosmic ray source models through secondary flux analysis

Abstract

We introduce the new version of SimProp, a Monte Carlo code for simulating the propagation of ultra-high energy cosmic rays in intergalactic space. This version, SimProp v2r4, together with an overall improvement of the code capabilities with a substantial reduction in the computation time, also computes secondary cosmogenic particles such as electron-positron pairs and gamma rays produced during the propagation of ultra-high energy cosmic rays. As recently pointed out by several authors, the flux of this secondary radiation and its products, within reach of the current observatories, provides useful information about models of ultra-high energy cosmic ray sources which would be hard to discriminate otherwise.

Figures

Tables

Table 1. Table 1 : Parameters of PSB cross sections [ 22 ] with Stecker and Salamon thresholds [ 23 ] (all energies in MeV)

		one-nucleon				two-nucleon
$A$	$Z$	${ϵ}_1^{\min }$	${ϵ}_1^0$	${\xi }_1$	${\Delta }_1$	${ϵ}_2^{\min }$	${ϵ}_2^0$	${\xi }_2$	${\Delta }_2$	$\zeta$
$56$	$26$	$10.2$	$18$	$0.98$	$8$	$18.3$	$22$	$0.15$	$7$	$0.95$
$55$	$25$	$8.1$	$18$	$0.93$	$7$	$17.8$	$24$	$0.20$	$8$	$0.95$
$54$	$26$	$8.9$	$18$	$0.93$	$7$	$15.4$	$24$	$0.20$	$8$	$0.95$
$53$	$24$	$7.9$	$18$	$1.03$	$7$	$18.4$	$24$	$0.10$	$8$	$0.95$
$52$	$24$	$10.5$	$18$	$1.08$	$7$	$18.6$	$24$	$0.05$	$8$	$0.95$
$51$	$23$	$8.1$	$19$	$1.02$	$7$	$19.0$	$25$	$0.11$	$6$	$0.95$
$50$	$24$	$9.6$	$19$	$1.03$	$8$	$16.3$	$25$	$0.10$	$6$	$0.95$
$49$	$22$	$8.1$	$19$	$1.03$	$8$	$19.6$	$25$	$0.10$	$6$	$0.95$
$48$	$22$	$11.4$	$19$	$1.03$	$8$	$19.9$	$25$	$0.10$	$6$	$0.95$
$47$	$22$	$8.9$	$19$	$1.03$	$8$	$18.7$	$25$	$0.10$	$6$	$0.95$
$46$	$22$	$10.3$	$19$	$1.03$	$8$	$17.2$	$25$	$0.10$	$6$	$0.95$
$45$	$21$	$6.9$	$19$	$0.97$	$9$	$18.0$	$26$	$0.15$	$8$	$0.95$
$44$	$20$	$11.1$	$20$	$0.92$	$9$	$19.1$	$26$	$0.20$	$8$	$0.96$
$43$	$20$	$7.9$	$20$	$0.97$	$8$	$18.2$	$26$	$0.15$	$8$	$0.96$
$42$	$20$	$10.3$	$20$	$1.02$	$7$	$18.1$	$26$	$0.10$	$8$	$0.96$
$41$	$19$	$7.8$	$20$	$0.92$	$6$	$17.7$	$26$	$0.20$	$8$	$0.96$
$40$	$20$	$8.3$	$20$	$0.84$	$6$	$14.7$	$26$	$0.28$	$10$	$0.96$
$39$	$19$	$6.4$	$20$	$0.73$	$7$	$16.6$	$25$	$0.38$	$12$	$0.98$
$38$	$18$	$10.2$	$18$	$0.86$	$8$	$18.6$	$22$	$0.24$	$8$	$0.98$
$37$	$17$	$8.4$	$20$	$0.81$	$7$	$18.3$	$24$	$0.28$	$7$	$1.00$
$36$	$18$	$8.5$	$22$	$0.82$	$12$	$14.9$	$22$	$0.25$	$12$	$1.00$
$35$	$17$	$6.4$	$20$	$0.87$	$7$	$17.3$	$26$	$0.22$	$10$	$1.00$
$34$	$16$	$10.9$	$22$	$0.87$	$12$	$20.1$	$22$	$0.20$	$12$	$1.00$
$33$	$16$	$8.6$	$22$	$0.82$	$12$	$17.5$	$22$	$0.25$	$12$	$1.00$
$32$	$16$	$8.9$	$22$	$0.97$	$12$	$16.2$	$30$	$0.10$	$12$	$1.00$
$31$	$15$	$7.3$	$21$	$0.85$	$8$	$17.9$	$29$	$0.20$	$12$	$1.02$

Table 2. Table 2 : Branching ratios for the number of nucleons ejected in the PSB model for  30 ​ MeV < ϵ ′ < 150 ​ MeV 30 MeV superscript italic-ϵ ′ 150 MeV 30~{}\mathrm{MeV}<\epsilon^{\prime}<150~{}\mathrm{MeV} as a function of the parent nucleus mass number

$A$	1 nuc.	2 nuc.	3 nuc.	etc.
3–4	$80\%$	$20\%$
9	$100\%$
10–22	$10\%$	$30\%$	$10\%$	$10\%$	$20\%$	$20\%$
23–56	$10\%$	$35\%$	$10\%$	$5\%$	$15\%$	$4.5\%$	$4.0\%$	$3.5\%$	$3.0\%$	$2.5\%$	$2.0\%$	$1.8\%$	$1.5\%$	$1.2\%$	$1.0\%$

Table 3. Table 3 : Source emissivities at z = 0 𝑧 0 z=0 used in figure 3

source	proton-only	mixed-composition scenario
evolution	scenario	soft sources	hard sources
none	13.8	19.2	0.43
SFR	06.2	08.7
AGN	03.4	04.7
all emissivities in ${10}^{45}\mathrm{erg}/{\mathrm{Mpc}}^3/\mathrm{yr}$

option	description	default
-h	prints the help and exits	none
-s	seed of the random number generator	65539
-N	number of events to be generated	100
-A	mass number of primary nuclei, $A_{\text{inj}}$ (chosen at random for each event with -A 0)	56
-e	${\log }_{10}\left(E_{\min }/\mathrm{eV}\right)$, where $E_{\min }$ is the minimum injection energy	17
-E	${\log }_{10}\left(E_{\max }/\mathrm{eV}\right)$, where $E_{\max }$ is the maximum injection energy	21
-g	injection spectral index, $\gamma$	1
-z	minimum source redshift, $z_{\min }$	0
-Z	maximum source redshift, $z_{\max }$	1
-r	distance between sources, $L_{\text{s}}$,666If $L_{\text{s}}$ is nonzero, then whenever an event is generated with it is rounded up to the next higher integer multiple of $\Delta z=H_0{L}_{\text{s}}=L_{\text{s}}/\left(4285\mathrm{Mpc}\right)$. in $\mathrm{Mpc}$	0
-L	EBL model, see section 4.1	1
-M	photodisintegration model, see section 4.2	0
-n	nucleon ejection scaling factor (only with -M 3 and -M 4)	1
-a	alpha-particle ejection scaling factor (only with -M 3 and -M 4)	1
-D	beta decay: 0 disabled, all nuclei treated as their respective beta-decay stable isobars; 1 enabled, treated as instantaneous	1
-S	treatment of pion production: -1 continuous energy loss approximation for protons, neglected for other nuclei (as in SimProp v2r0); 0 continuous energy loss for both protons and other nuclei; 1 stochastic, on the CMB only; 2 stochastic, on both the CMB and the EBL	1
-p	electrons and positrons: 0 disregarded, 1 individually written to output file (warning: results in very large output files), 2 binned according to production redshift (by default $[{10}^{-4.0},{10}^{-3.8}),\dots ,[{10}^{+0.8},{10}^{+1.0})$, but can be changed in the file src/Output.h before compiling), 3 only total energy written	0
-o	output type (see below): 0 old (nuc and ev trees); 1 new (summary tree); 2 both	0

branch	description
evt	event number (starting from 0)
branch	branch generation number: $0$ for the primary; incremented by $1$ from the parent branch in stochastic interactions and decays
intmult	$0$ if the particle reaches Earth; $1$ if the particle undergoes pion production or disintegration with one-nucleon ejection; $n$ if it undergoes $n$-nucleon ejection; $4$ if it undergoes alpha-particle ejection; $1000$ for photons and electrons, whose propagation is not implemented in SimProp and can be simulated with external tools; $1000+n$ if the particle decays into $n$ particles
Acurr	mass number during the branch (0 for leptons, photons and pions)
Zecurr	electric charge during the branch
Flav	flavours of leptons ($+1$ for $e^{-}$ and ${\nu }_e$, $-1$ for $e^{+}$ and ${\overline{\nu }}_e$, $+2$ for ${\mu }^{-}$ and ${\nu }_{\mu }$, $-2$ for ${\mu }^{+}$ and ${\overline{\nu }}_{\mu }$) at production, $0$ for all other particles
zOri	redshift at the beginning of the branch
zEnd	redshift at the end of the branch
EOri	energy at the beginning of the branch, in $\mathrm{eV}$
EEnd	energy at the end of the branch, in $\mathrm{eV}$
Dist	comoving distance travelled (see appendix B), in $\mathrm{Mpc}$

branch	description
timexev	CPU time used during the event, in seconds
branxev	total number of branches in the event
seed	seed of the random number generator at the end of the event

branch	description
event	event number (starting from 0)
injEnergy	injection energy, in $\mathrm{eV}$
injRedshift	source redshift
injDist	source comoving distance (see appendix B), in $\mathrm{Mpc}$
injA	injection mass number
injZ	injection atomic number
nNuc	number of nuclei reaching Earth
nucEnergy[nNuc]	energies of nuclei reaching Earth, in $\mathrm{eV}$
nucA[nNuc]	mass numbers of nuclei reaching Earth
nucZ[nNuc]	atomic numbers of nuclei reaching Earth
nPho	number of photons produced
phoEProd[nPho]	production energies of photons, in $\mathrm{eV}$
phozProd[nPho]	redshifts of production points of photons
nEle	when using -p 0: 0 (electrons and positrons disregarded); when using -p 1: number of electrons and positrons produced; when using -p 2: number of $e^{\pm }$ production redshift bins used; when using -p 3: 1 (only the totals are given)
eleZ[nEle]	$-1$ for electrons, $+1$ for positrons (with -p 2: sums over $z_{\text{prod}}$ bins; with -p 3: sum over all $e^{\pm }$)
eleEProd[nEle]	production energies of $e^{\pm }$, in $\mathrm{eV}$ (with -p 2: sums over $z_{\text{prod}}$ bins; with -p 3: sum of $E_{\text{prod}}/\left(1+z_{\text{prod}}\right)$ over all $e^{\pm }$)
elezProd[nEle]	redshifts of production points of $e^{\pm }$ (with -p 2: logarithmic centres of $z_{\text{prod}}$ bins; with -p 3: 0)
nNeu	number of neutrinos reaching Earth
neuEnergy[nNeu]	energies of neutrinos reaching Earth, in $\mathrm{eV}$
neuFlav[nNeu]	flavours of neutrinos at production ($+1$ for ${\nu }_e$, $-1$ for ${\overline{\nu }}_e$, $+2$ for ${\nu }_{\mu }$, $-2$ for ${\overline{\nu }}_{\mu }$)

Table 8. Table 8 : Computation times required by SimProp  v2r4, using the PSB photodisintegration model ( -M 0 ) and the Gilmore et al. 2012 fiducial EBL model ( -L 7 ). Similar results are obtained with the TALYS photodisintegration model ( -M 4 < xsect_Gauss2_TALYS-restored.txt ) and/or the Domínguez et al. 2011 best-fit EBL model ( -L 4 ).

$z_{\text{inj}}$	${\log }_{10}\left(E_{\text{inj}}/\mathrm{eV}\right)\to$	$[18,19)$			$[19,20)$
$\downarrow$	$A_{\text{inj}}\to$	$1$	$14$	$56$	$1$	$14$	$56$
$[0.0,0.2)$	-S 0 -p 0	0.06	5.92	13.70	0.05	6.29	18.08
	-S 1 -p 0	1.13	9.88	23.13	1.05	12.12	30.84
	-S 2 -p 0	3.24	21.78	50.62	2.93	27.05	63.37
	-S 2 -p 2	13.78	126.90	444.78	13.69	181.69	574.90
$[0.2,0.5)$	-S 0 -p 0	0.09	9.83	21.15	0.08	8.62	26.40
	-S 1 -p 0	2.41	22.70	53.72	2.27	26.24	70.07
	-S 2 -p 0	6.49	48.32	129.47	5.55	63.66	173.37
	-S 2 -p 2	57.81	633.57	1134.89	55.74	873.60	1775.10
$[0.5,1.0)$	-S 0 -p 0	0.10	16.95	29.82	0.13	10.33	32.94
	-S 1 -p 0	3.25	47.39	99.87	3.89	46.48	134.69
	-S 2 -p 0	9.36	110.51	254.83	9.87	117.60	341.75
	-S 2 -p 2	97.62	1407.66	2269.57	126.48	2109.06	3936.28
all times in milliseconds per event

Equations

$$n_{\text{pairs}}$$

$$\frac{1}{\tau }=\frac{1}{2{\Gamma }^2}{\int }_{{ϵ}_{\min }^{\prime }}^{2\Gamma ϵ}{\int }_{ϵ=0}^{+\infty }\frac{n_{\gamma }(ϵ)}{{ϵ}^2}\text{differential}ϵ\sigma ({ϵ}^{\prime }){ϵ}^{\prime }\text{differential}{ϵ}^{\prime },$$

$$I(ϵ)$$

$$\frac{1}{\tau }=\frac{1}{4m^2{\Gamma }^2}{\int }_{{ϵ}_{\min }^{\prime }/2\Gamma }^{+\infty }\Phi (m^2+4m\Gamma ϵ)\frac{n_{\gamma }(ϵ)}{2{ϵ}^2}\text{differential}ϵ=\frac{1}{{\Gamma }^2}{\int }_{{ϵ}_{\min }^{\prime }}^{+\infty }I\left(\frac{{ϵ}^{\prime }}{2\Gamma }\right){ϵ}^{\prime }\sigma ({ϵ}^{\prime })\text{differential}{ϵ}^{\prime }.$$

$$n_{\text{CMB}}(ϵ)$$

$${\tau }_{\text{pion,CMB}}^{-1}(A,\Gamma ,z)=(1+z{)}^{3}A{\tau }_{\text{pion,CMB}}^{-1}(\text{proton},(1+z)\Gamma ,z=0),$$

$$p(ϵ)\text{differential}ϵ$$

$$p(s|ϵ)\text{differential}s$$

$$E_{\pi }^{\ast }$$

$$\cos {\theta }_{12}=\frac{p_{\mathrm{e}}^{\ast 2}-p_{{\nu }_1}^{\ast 2}-p_{{\nu }_2}^{\ast 2}}{2p_{{\nu }_1}^{\ast }{p}_{{\nu }_2}^{\ast }};$$

$${\sigma }_i({ϵ}^{\prime })=\{\begin{array}{ll}{\xi }_i\frac{{\Sigma }_{\text{d}}}{W_i}\exp (-\frac{({ϵ}^{\prime }-{ϵ}_i^0{)}^2}{{\Delta }_i^2}),& {ϵ}_i^{\min }<{ϵ}^{\prime }<{ϵ}_1,i=1,2;\\ b_i\frac{\zeta {\Sigma }_{\text{d}}}{{ϵ}_{\max }-{ϵ}_1},& {ϵ}_1<{ϵ}^{\prime }<{ϵ}_{\max },i=1,2,\dots ;\end{array}$$

$$W_i$$

$${\sigma }_i({ϵ}^{\prime })=\{\begin{array}{ll}{\xi }_i\text{quantity}(1+\frac{({ϵ}^{\prime }-{ϵ}_i^0{)}^2}{{\Delta }_i^2}{)}^{-1},& {ϵ}_i^{\min }<{ϵ}^{\prime }<{ϵ}_1,i=1,2;\\ \text{same as for models 0 and 1},& {ϵ}_1<{ϵ}^{\prime }<{ϵ}_{\max }.\end{array}$$

$${\sigma }_i({ϵ}^{\prime })=\{\begin{array}{ll}{h}_i^1\text{quantity}(1+\frac{({ϵ}^{\prime }-x_i^1{)}^2}{{w_i^1}^2}{)}^{-1}+h_i^2\text{quantity}(1+\frac{({ϵ}^{\prime }-x_i^2{)}^2}{{w_i^2}^2}{)}^{-1},& t_i<{ϵ}^{\prime }<{ϵ}_1,i=N,\alpha ;\\ c_i,& {ϵ}_1<{ϵ}^{\prime }<{ϵ}_{\max },i=N,\alpha .\end{array}$$

$${\sigma }_i({ϵ}^{\prime })=\{\begin{array}{ll}{h}_i^1\exp (-\frac{({ϵ}^{\prime }-x_i^1{)}^2}{w_i^1}),& t_i<{ϵ}^{\prime }<{ϵ}_1,i=N,\alpha ;\\ c_i,& {ϵ}_1<{ϵ}^{\prime }<{ϵ}_{\max },i=N,\alpha .\end{array}$$

$${\sigma }_N$$

$${\sigma }_{\alpha }$$

$$d_{\text{C}}(z)={\int }_0^z\frac{\text{differential}z}{H_0\sqrt{(1+z{)}^3{\Omega }_{\text{m}}+{\Omega }_{\Lambda }}}.$$

$$d_{\text{T}}(z)={\int }_0^z\text{absolutevalue}\text{derivative}tz\text{differential}z={\int }_0^z\frac{\text{differential}z}{H_0(1+z)\sqrt{(1+z{)}^3{\Omega }_{\text{m}}+{\Omega }_{\Lambda }}}.$$

$$L={\int }_0^{+\infty }EQ(E_{\text{inj}})\text{differential}{E}_{\text{inj}},$$

$$J_i(E_{\text{Earth}})=\frac{c}{4\pi }{\int }_{t_{\min }}^{t_0}{\int }_0^{+\infty }{T}_{ij}(E_{\text{Earth}}|E_{\text{inj}},t){\mathcal{Q}}_j(E_{\text{inj}},t)\text{differential}{E}_{\text{inj}}\text{differential}t,$$

Full text

1 Introduction

Cosmic rays (high-energy charged particles from outer space, mainly protons and other nuclei) with energies up to a few hundred EeV ($1~{}\mathrm{EeV}=10^{18}~{}\mathrm{eV}\approx 0.16~{}\mathrm{J}$) have been observed since the 1960s. Half a century later, the origin of such particles is still uncertain, but there is now a wide consensus that most of the cosmic rays with energies above a few EeV (known as ultra-high-energy cosmic rays, UHECRs) originate from extragalactic sources [1]. If this is the case, the energy, type and trajectory of propagating particles can change in their intergalactic journey through different processes. These include: (i) the adiabatic energy loss that all particles travelling cosmological distances undergo due to the expansion of the Universe (the redshift loss); (ii) photonuclear interactions with cosmic microwave background (CMB) and infrared/visible/ultraviolet extragalactic background light (EBL) photons; and (iii) deflections by intergalactic and galactic magnetic fields. Photonuclear interactions include electron–positron pair production, disintegration of nuclei (whereby one or more nucleons or other light fragments are ejected from a nucleus), and photohadronic interactions producing one or more mesons (mainly pions). These processes need to be taken into account when predicting the energy spectrum, mass composition and distribution of arrival directions of UHECRs expected at the Earth assuming a given model for their sources. In addition, photonuclear processes may result in the production of secondary particles such as nucleons, electron–positron pairs, neutrinos and gamma rays. The detection of this radiation or lack thereof can further help in constraining scenarios for the sources of UHECRs.

Some of the quantities relevant to UHECR propagation are known with good precision. These include: (i) the redshift energy loss rate (which in the standard $\Lambda$CDM cosmology can be computed from the FLRW metric as $-\quantity(\differential{\ln E}/\differential{t})_{\text{ad}}=H(z)=H_{0}\sqrt{(1+z)^{3}\Omega_{\text{m}}+\Omega_{\Lambda}}$, where $H_{0}\approx 70~{}\mathrm{km}/\mathrm{s}/\mathrm{Mpc}$, $\Omega_{\text{m}}\approx 0.3$ and $\Omega_{\Lambda}\approx 0.7$), (ii) the CMB spectrum and its cosmological evolution (a blackbody spectrum with temperature $T=(1+z)T_{0}$, $T_{0}\approx 2.7~{}\mathrm{K}$), (iii) photonuclear cross sections for pair production (computed from the Bethe–Heitler formula) and pion production (measured in accelerator-based experiments). Other quantities are still considerably uncertain, and phenomenological models are needed to estimate them; these include: (i) the spectrum and cosmological evolution of the EBL (especially in the far infrared and at high redshifts) and (ii) exclusive cross sections for certain photodisintegration channels (especially those where charged fragments are ejected). These uncertainties have a non-negligible impact on UHECR propagation studies as widely discussed in ref. [2].

SimProp is a simple and fast Monte Carlo code for simulating UHECR propagation through intergalactic space. The processes taken into account are the redshift energy loss, photonuclear interactions with background photons, and decays of unstable particles (pions, muons and unstable nuclei). The photon backgrounds taken into account in SimProp are the CMB and the EBL, with several models available for the EBL spectrum and evolution (section 4.1). The photonuclear interaction processes implemented in SimProp are the electron–positron pair production, approximated as a continuous energy loss process; the disintegration of nuclei, for which different models are available (section 4.2); and photohadronic interactions, all approximated as single-pion production with branching ratios computed assuming isospin invariance, distribution of outgoing pion directions taken to be isotropic in the centre-of-mass (CoM) frame, and nuclei treated as collections of free nucleons. Particle decays are approximated as instantaneous, as decay lengths are generally much shorter than all other relevant length scales. Intergalactic magnetic fields are neglected in SimProp; consequently, UHECR propagation through intergalactic space is taken to be rectilinear.

The main new feature of SimProp v2r4 compared to the previous versions [3, 4, 5, 6] is the computation of the production rates of high-energy secondary electrons, positrons and gamma rays produced during the propagation of UHECRs. Such particles can initiate electromagnetic cascades via inverse Compton scattering and electron–positron pair production against photons of the CMB, EBL and radio backgrounds, resulting in a diffuse background of photons up to TeV energies. Measurements of the diffuse extragalactic gamma-ray background can provide an upper bound to the energy density of such cascades, which can be used to constrain models of UHECR sources [7, 8, 9, 10].

This paper is structured as follows: in section 2 we briefly recall the history of SimProp versions; in section 3 we describe technical details of the implementation of SimProp v2r4; in section 4 we list and discuss the various models of extragalactic background light and photodisintegration cross sections available; in section 5 we discuss the new capabilities of SimProp v2r4 and how they can be useful for testing UHECR source models. In appendix A we describe the command-line options SimProp v2r4 accepts and the structure of its output files. In appendix B we define the various measures of cosmological distances used in the text. Finally, in appendix C we report the computation time required by SimProp simulations in a few selected situations.

Unless otherwise specified, natural units $c=\hbar=1$ are used.

SimProp is available upon request to SimProp-dev@aquila.infn.it.

2 Motivation and history

SimProp was originally developed as a refinement of the analytic models by Aloisio, Berezinsky and Grigorieva [11, 12] in order to have a publicly available Monte Carlo code for the community to use, at a time when most UHECR propagation studies used closed-source simulation codes such as that by Allard et al. [13] or Hooper et al. [14].

Since the first SimProp release [3], more public codes have become available, including CRPropa [15, 16], Hermes [17], and TransportCR [18]. The usefulness of having several independently developed codes is that it allows us to cross-check results, improving the computation reliability. For example, through the comparison of the results obtained with CRPropa and SimProp it was possible to assess the relevance of uncertainties in a few astrophysical and nuclear physics quantities [2], neglected before. The same comparison brought to the conclusion that CRPropa and SimProp show a substantial agreement [2], confirming the validity of these computation schemes.

Another test of the correctness of SimProp was the comparison of Engel et al. [19] neutrino fluxes to SimProp results in the same physical scenario222Note that the model considered there assumes a pure proton injection with practically no injection cutoff, which is no longer considered plausible due to UHECR composition data from the Pierre Auger Observatory [20], and neglects interactions with the EBL, which are the dominant origin of cosmogenic neutrinos below about 0.1 EeV [2]. showing very good agreement except in the electron antineutrino flux at the highest energies (which is in any case subdominant) due to the approximation of all photohadronic interactions as single-pion production in SimProp [5]. Finally, a new systematic comparison of the outputs of SimProp and CRPropa mainly in terms of secondary neutrinos and gamma rays is planned and will be soon released [21].

In the first released version, SimProp v2r0 [3], the only processes treated stochastically were the sampling of the injection redshift, the energy of primary particles and the photodisintegration of nuclei. All other processes, namely the adiabatic energy loss, pair production and pion production, were treated deterministically, as in the analytic models SimProp was based on [11, 12]. Pion production was only taken into account for protons interacting with CMB photons and approximated as a continuous energy loss. Photodisintegration was treated according to the Puget–Stecker–Bredekamp (PSB) model [22] as refined by Stecker and Salamon [23], taking into account the EBL using the Stecker et al. model [24] or a power-law approximation thereof as well as the CMB.

Starting from the following version, SimProp v2r1 [4], the pion production process is also treated stochastically, for both free protons and nucleons bound within nuclei, in order to compute fluxes of $\mathrm{EeV}$ secondary neutrinos thereby produced, also assessing the accuracy of the continuous energy loss approximation for that process [2].

SimProp v2r2 [5], apart from fixing a few bugs in SimProp v2r1, added an option to also take into account pion production on the EBL, in order to compute secondary neutrino fluxes down to $\mathrm{PeV}$ energies [25], also including the choice of using the Kneiske et al. model [26] for the spectrum and cosmological evolution of the EBL.

Starting from SimProp v2r3 [6, 27], the user can also choose the EBL models from Domínguez et al. [28] (best-fit, lower-limit or upper-limit) or Gilmore et al. [29] (fiducial), and one of four different parametrizations for photodisintegration cross sections with user-defined parameter values. Processes where alpha particles are ejected (e.g. ${{}^{A}Z}+\gamma\to{{}^{A-4}(Z-2)}+{{}^{4}\mathrm{He}}$) are also implemented, and their cross sections can be scaled by a user-defined factor via a simple command-line option, allowing the user to assess the effects of these poorly known quantities on the results.

In the present version, SimProp v2r4, secondary electron–positron pairs and photons produced by the propagation of UHECRs can also be written in the output file, so that external tools, such as DINT [30], ELMAG [31], or $\mathit{EleCa}$ [32] can then be used to simulate the electromagnetic cascades initiated by such particles interacting with CMB, EBL and radio backgrounds. In addition, the user can now specify an arbitrary spectral index for the energy distribution of the injected particles. Finally, a few adjustments and bug fixes have been applied to the new version of SimProp, reducing the overall computation time; a few examples of the computation times required by the current SimProp version are reported in appendix C.

3 Structure of the code

SimProp is written in C++, and makes use of features from the ROOT [33] framework. The random number generator used in SimProp is TRandom3 from ROOT, which is based on the Mersenne twister [34].

A SimProp run consists of $N$ events, each including the generation of a primary particle with mass number $A_{\text{inj}}$, initial energy $E_{\text{inj}}$ sampled from a power-law distribution from $E_{\min}$ to $E_{\max}$ with spectral index $\gamma$ (i.e. $\differential{N}\propto E_{\text{inj}}^{-\gamma}\differential{E_{\text{inj}}}$), and source redshift $z_{\text{inj}}$ uniformly distributed from $z_{\min}$ to $z_{\max}$.333This can be converted to different distributions by weighing each event by an appropriate function of $z_{\text{inj}}$, for example by $w(z_{\text{inj}})\propto{1}/{\sqrt{(1+z_{\text{inj}})^{3}\Omega_{\text{m}}+\Omega_{\Lambda}}}$ for a constant density of sources per unit comoving volume (see appendix B). The propagation of particles is followed, along with that of any secondary particles produced during propagation. The injection parameters $N$, $A_{\text{inj}}$, $E_{\min}$, $E_{\max}$, $\gamma$, $z_{\min}$ and $z_{\max}$ can be specified via command-line options (appendix A).

The interval between the production of a particle (at injection for primaries, or at an interaction for secondaries) and its decay, stochastic interaction or arrival at the Earth is called a branch. During a branch, a particle keeps the same mass number and electric charge, but loses energy through continuous processes, such as adiabatic and pair production energy loss. Each of the outgoing particles in a stochastic interaction or decay starts a new branch, even if it is of the same type as the parent (e.g. as in neutral pion production). Secondary electrons and positrons from pair production also start a new branch, but the parent nucleus continues the same branch.

3.1 Propagation of protons and stable nuclei

When propagating a proton or a stable nucleus, the redshift interval between its production point $z_{\text{prod}}$ and [math] (Earth) is divided into steps $z_{0}=z_{\text{prod}}$, $z_{1}$, …, $z_{n}=0$. During each step, there are two types of processes the particle can undergo:

• •

those treated in the approximation of continuous (deterministic) energy losses, namely the adiabatic energy loss, pair production on CMB photons, and (if the option -S -1 or -S 0 is used) pion production on CMB photons;

• •

those treated as discrete (stochastic) interactions, with the interaction point and the energy, type and/or number of outgoing particles to be sampled stochastically, namely the photodisintegration of nuclei and pion production (on CMB photons by default, if the option -S 1 is used, or also on EBL photons if -S 2 is used).

3.1.1 Continuous energy losses

At each redshift step $(z_{i-1},z_{i}]$, the processes computed in the continuous energy loss approximation are treated by numerically integrating the differential equation for the logarithm of the Lorentz factor $\ln\Gamma$ from $z_{i-1}$ to $z_{i}$, as described in ref. [3]. The values used for the cosmological parameters are $H_{0}=7.17\times 10^{-11}~{}\mathrm{yr}^{-1}\approx 70~{}\mathrm{km}/\mathrm{s}/\mathrm{Mpc}$, $\Omega_{\text{m}}=0.3$ and $\Omega_{\Lambda}=0.7$. The pair-production energy loss rates for protons at $z=0$ are interpolated from a table computed as described in ref. [35]. For nuclei heavier than protons, pair production energy losses are approximated as those of protons multiplied by the factor $A/4$ (which is typically equal or very close to the exact factor $Z^{2}/A$).

The number of secondary electron–positron pairs generated in the step and optionally written to the output file is sampled from a Poisson distribution with mean

[TABLE]

on the assumption that each pair is generated with energy $K$ times that of the parent particle, as it is at the kinematic threshold. The production points of the secondary pairs are chosen at random uniformly from $z_{i-1}$ to $z_{i}$, and the energies of the electrons and positrons as $KAm_{N}\Gamma/2$, where $\log\Gamma$ is chosen at random uniformly from $\log\Gamma_{i-1}$ to $\log\Gamma_{i}$.

The approximation that the fraction of energy carried by each pair is $K$ is not very accurate for energies much higher than the threshold. Nevertheless, high-energy electrons produced at even moderate $z$ result in cascading photons at TeV energies. This radiation has a spectrum quite independent of the individual energies of the initial electrons [36]; the only quantity with observable consequences is the total energy emitted in pairs at each interval, which is not affected by the approximation above.

3.1.2 Discrete interactions

An extension of the scheme described in ref. [3] is used in order to determine whether a particle undergoes a discrete interaction within a given redshift step and, if it does, to sample the interaction point within the step, type, number and energies of the interaction products.

The total interaction rate (probability per unit time) $1/\tau$ of an ultrarelativistic particle with background photons is given by

[TABLE]

where $\epsilon^{\prime}_{\min}$ is the threshold energy of the process at hand, $n_{\gamma}(\epsilon)$ is the density of background photons, $\sigma(\epsilon^{\prime})$ is the interaction cross section, $\Gamma$ is the particle’s Lorentz factor and all primed quantities are intended in the CoM reference frame. Since $1/\tau$ is linear in both the background photon density and the cross sections, it can be computed as the sum of the interaction rates for each type of process and photon background, $\tau^{-1}=\sum_{ij}\tau_{ij}^{-1}$. When the particle is found to have interacted, the probability that the interaction was of type $i$ with a photon of the $j$-th background is then $p_{ij}=\tau_{ij}^{-1}/\tau^{-1}$.

If we introduce the quantities

[TABLE]

where $m$ is the mass of the particle (a nucleus in the case of disintegration and a nucleon in the case of pion production) and $s$ is the CoM energy squared $s=m^{2}+2m\epsilon^{\prime}$, then eq. (3.2) can be rewritten as

[TABLE]

The photon background $n_{\gamma}$ is the sum of two terms, one for the CMB and one for the EBL. The CMB spectrum is precisely known at all redshifts, being a blackbody spectrum:

[TABLE]

where $T=(1+z)T_{0}$, and the value of $k_{\text{B}}T_{0}$ used in SimProp is $0.2327~{}\mathrm{meV}$ ($T_{0}\approx 2.7~{}\mathrm{K}$). This implies that the $z$-dependence of the interaction lengths on the CMB is given by $\tau_{\text{CMB}}^{-1}(\Gamma,z)=(1+z)^{3}\tau^{-1}_{\text{CMB}}((1+z)\Gamma,z=0)$. On the other hand, EBL is not precisely known, and needs to be approximated using phenomenological models; those available in SimProp are listed in section 4.1.

As for the cross sections, we used the total photohadronic cross section computed by SOPHIA [37] for protons as the pion photoproduction cross section $\sigma_{\text{pion}}(\epsilon^{\prime})$, and numerically integrated the second integral in eq. (3.3) to obtain a table of values from which we interpolate $\Phi_{\text{pion}}(s)$. As for photodisintegration, various models are available, listed in section 4.2.

From these quantities, $\tau_{\text{pion,CMB}}^{-1}(A,\Gamma,z)$ is computed as

[TABLE]

where $\tau^{-1}_{\text{pion,CMB}}$ for protons at $z=0$ is interpolated from a table with values obtained by numerically integrating the first integral in eq. (3.4). $\tau^{-1}_{\text{pion,EBL}}$ is computed as a function of $\Gamma$ and $z$ via a 2D interpolation from a table obtained by numerically integrating the second integral in eq. (3.4). As for pion production, we assume that a nucleus behaves as $A$ independent nucleons, i.e., $\tau_{\text{pion}}^{-1}(A,\Gamma,z)=A\tau_{\text{pion}}^{-1}(\text{proton},\Gamma,z)$ because the energies involved are much larger than the binding energy per nucleon. Finally, if needed, $\tau^{-1}_{\text{disi}}$ is computed by numerically integrating the second integral in eq. (3.4). Then we have $\tau^{-1}=\tau_{\text{pion,CMB}}^{-1}+\tau_{\text{pion,EBL}}^{-1}+\tau_{\text{disi}}^{-1}$.

After a particle is found to have interacted and the type of interaction is determined, the number, type and energy of the outgoing particles are sampled as follows.

Photodisintegration.

When a nucleus is photodisintegrated, one of the channels implemented in the photodisintegration model used (see section 4.2) is chosen at random, with probabilities proportional to their respective interaction rates. The energy of the parent nucleus is assumed to be split among the residual nucleus and the ejected fragments proportionally to their mass, i.e. all the fragments inherit the Lorentz factor of the original nucleus.

Pion photoproduction.

When a pion is photoproduced, if the incoming particle is a nucleus, the nucleon that undergoes the interaction is chosen at random. We approximate all photohadronic processes as single-pion production; assuming isospin invariance, a neutral pion is produced ($\mathrm{p}+\gamma\to\mathrm{p}+\pi^{0}$, $\mathrm{n}+\gamma\to\mathrm{n}+\pi^{0}$) with probability $2/3$ and a charged pion is produced ($\mathrm{p}+\gamma\to\mathrm{n}+\pi^{+}$, $\mathrm{n}+\gamma\to\mathrm{p}+\pi^{-}$) with probability $1/3$.

In order to sample the pion energy, first the photon energy $\epsilon$ in the lab frame is sampled from its marginal distribution444In practice, we use the fact that the marginal distribution of $\epsilon$ corresponds to a distribution of $I(\epsilon)$ proportional to $\Phi(m^{2}+4m\Gamma\epsilon)$, and the conditional distribution of $s$ given $\epsilon$ corresponds to a uniform distribution of $\Phi(s)$, so we actually sample $I$ and $\Phi$ and invert the functions to find the corresponding $\epsilon$ and $s$.

[TABLE]

where $m$ is the nucleon mass and $\epsilon^{\prime}_{\min}=m_{\pi}+m_{\pi}^{2}/2m$; then, given $\epsilon$, the squared CoM energy $s$ is sampled from its conditional distribution

[TABLE]

from $s$ the pion energy and momentum in the CoM frame are calculated as

[TABLE]

and the Lorentz factor of the transformation from the CoM frame to the lab frame as $\gamma=m\Gamma/\sqrt{s}$; then the pion energy is converted to the lab frame as $E_{\pi}=\gamma(E^{*}_{\pi}+p^{*}_{\pi}\cos\theta_{\pi})$, where the distribution of $\theta_{\pi}$ is approximated as isotropic ($\cos\theta_{\pi}$ uniformly distributed from $-1$ to $1$). In the lab frame, the momentum component orthogonal to the original direction is much smaller than that parallel to it, by a factor of order $\epsilon/E\sim 10^{-20}$, so we neglect transverse components assuming one dimensional propagation.

The outgoing particles are then the pion with energy $E_{\pi}$, the nucleon with energy $m\Gamma-E_{\pi}$, and (in the case of nuclei) the residual nucleus with mass number $A-1$ and energy $(A-1)m\Gamma$.

3.2 Decay of unstable particles

When an unstable particle is produced, it is assumed to decay instantaneously, as decay lengths are generally much shorter than all other relevant length scales [11, 12]. The energies of the decay products are sampled as follows.

Beta decay of neutrons and unstable nuclei.

Neutrons and nuclei not in the list of beta-decay stable isobars are assumed to immediately undergo beta decay. The $Q$-value of the reaction is read from a table taken from ref. [38] or, for nuclei not on that table, estimated via the semi-empirical mass formula.

Then, the electron energy in the nucleus rest frame $E^{*}_{\mathrm{e}}$ is sampled from a distribution $\propto(E^{*2}_{\mathrm{e}}-m_{\mathrm{e}}^{2})^{1/2}E^{*}_{\mathrm{e}}(Q-(E^{*}_{\mathrm{e}}-m_{\mathrm{e}}))^{2}$ (i.e., neglecting electromagnetic effects) and the neutrino energy is calculated as $E^{*}_{\nu}=Q-(E^{*}_{\mathrm{e}}-m_{\mathrm{e}})$; the recoil of the nucleus is neglected. The electron and neutrino energies are converted to the lab frame by $E_{\mathrm{e}}=\Gamma E^{*}_{\mathrm{e}}(1-\cos\theta_{\mathrm{e}})$ and $E_{\nu}=\Gamma E^{*}_{\nu}(1-\cos\theta_{\nu})$, where $\Gamma$ is the Lorentz factor of the nucleus, and $(1-\cos\theta_{\mathrm{e}})$ and $(1-\cos\theta_{\nu})$ are sampled from the uniform distribution from [math] to $2$.

The decay products are then the daughter nucleus (with the same energy and mass number $A$ as the parent, with electric charge $Z$ incremented in $\beta^{-}$ decay and decremented in $\beta^{+}$ decay), the neutrino ($\bar{\nu}_{e}$ in $\beta^{-}$ decay, $\nu_{e}$ in $\beta^{+}$ decay) and the electron or positron.

Neutral pion decay.

A $\pi^{0}$ with energy $E_{\pi}$ decays into two photons with energy $E_{\gamma_{1}}$ distributed uniformly from 0 to $E_{\pi}$ and $E_{\gamma_{2}}=E_{\pi}-E_{\gamma_{1}}$.

Charged pion decay.

A $\pi^{\pm}$ with energy $E_{\pi}$ decays into a muon with energy $E_{\mu}$ distributed uniformly from [math] to $(1-m_{\mu}^{2}/m_{\pi}^{2})E_{\pi}$ and a neutrino with energy $E_{\nu}=E_{\pi}-E_{\mu}$.

Muon decay.

A muon with energy $E_{\mu}$ decays into two neutrinos and an electron; their energies are sampled as follows:

• •

the energies of the neutrinos in the muon rest frame $E^{*}_{\nu_{1}}$ and $E^{*}_{\nu_{2}}$ are sampled independently uniformly from 0 to $m_{\mu}/2-m_{\mathrm{e}}^{2}/2m_{\mu}$, and that of the electron is $E^{*}_{\mathrm{e}}=m_{\mu}-E^{*}_{\nu_{1}}-E^{*}_{\nu_{2}}$;

• •

the corresponding momenta are computed as $p^{*}_{\nu_{1}}=E^{*}_{\nu_{1}}$, $p^{*}_{\nu_{2}}=E^{*}_{\nu_{2}}$, $p^{*}_{\mathrm{e}}=\sqrt{E^{*2}_{\mathrm{e}}-m_{\mathrm{e}}^{2}}$;

• •

if these values violate any of the constraints $E^{*}_{\mathrm{e}}\geq m_{\mathrm{e}}$, $p^{*}_{\nu_{1}}\leq p^{*}_{\nu_{2}}+p^{*}_{\mathrm{e}}$, $p^{*}_{\nu_{2}}\leq p^{*}_{\mathrm{e}}+p^{*}_{\nu_{1}}$, or $p^{*}_{\mathrm{e}}\leq p^{*}_{\nu_{1}}+p^{*}_{\nu_{2}}$, they are discarded and a new $E^{*}_{\nu_{1}},E^{*}_{\nu_{2}}$ pair is sampled;

• •

the angle $\theta_{12}$ between the two neutrinos is given by

[TABLE]

• •

the angle $\theta_{1}$ between the first neutrino and the line of sight is isotropic, i.e. $\cos\theta_{1}$ uniform from $-1$ to $1$;

• •

the angle $\phi$ between the second neutrino and the plane containing the line of sight and the first neutrino is uniform from 0 to $2\pi$;

• •

the angle $\theta_{2}$ between the second neutrino and the line of sight is given by $\cos\theta_{2}=\cos\theta_{12}\cos\theta_{1}-\sin\theta_{12}\sin\theta_{1}\cos\phi$;

• •

finally, the neutrino energies are transformed to the lab frame via $E_{\nu_{1}}=\gamma(E^{*}_{\nu_{1}}+p^{*}_{\nu_{1}}\cos\theta_{1})$ and $E_{\nu_{2}}=\gamma(E^{*}_{\nu_{2}}+p^{*}_{\nu_{2}}\cos\theta_{2})$, and the electron energy is computed as $E_{\mathrm{e}}=E_{\mu}-E_{\nu_{1}}-E_{\nu_{2}}$.

3.3 Other particles: neutrinos, photons and electron–positron pairs

The propagation of neutrinos is trivial: no interaction is possible and the only energy loss is the redshift loss, so a neutrino produced with energy $E_{\text{prod}}$ at redshift $z_{\text{prod}}$ will reach Earth with energy $E_{\text{Earth}}=E_{\text{prod}}/(1+z_{\text{prod}})$. Neutrino flavours are recorded only at the production; flavour oscillations are not implemented in the code, as it is trivially easy to calculate flavour transformations from the source to the observer555For example, assuming tribimaximal mixing, a pure $\nu_{\mathrm{e}}$ flux at production results in $\nu_{\mathrm{e}}:\nu_{\mu}:\nu_{\tau}=5:2:2$ ratios at Earth, a pure $\nu_{\mu}$ or $\nu_{\tau}$ flux at production results in $4:7:7$ ratios at Earth, and $1:2:0$ ratios at production result in $1:1:1$ ratios at Earth..

The energy and production redshift of secondary photons or electron–positron pairs are recorded in the output files so that the electromagnetic cascades initiated by them can be simulated with external tools such as DINT [30], ELMAG [31], or $\mathit{EleCa}$ [32], as the propagation of photons and electron–positron pairs is not yet implemented in SimProp. Since the spectral shape of electromagnetic cascades at the Earth does not depend on the initial energy (or initiating particle type) and depends only weakly on the initial redshift [36], the expected flux of secondary gamma rays observed at the Earth can be estimated from the total energy deposited in photons and pairs in each redshift interval.

4 Models of physical quantities

4.1 Extragalactic background light

It is very difficult to directly measure the EBL because the foreground, the zodiacal light, is orders of magnitude larger than the EBL itself. Therefore, phenomenological models are used to describe its spectrum and its evolution with redshift, based on the time evolution of galaxy populations as estimated using different approaches. The results of these methods are not always compatible each other, even in the case of recent models, especially for long wavelengths (far IR) and at large redshifts.

The following EBL models are available in SimProp, and can be chosen via the -L command-line option:

[TABLE]

The corresponding spectral energy densities are plotted in figure 1. Note that EBL models up to -L 3 are included for backward compatibility, but they are in tension with Fermi gamma- ray observations of blazars [39], and their use is not recommended except for testing purposes.

In models 1, 4, 5, 6 and 7, $n_{\text{EBL}}(\epsilon,z)$ is found via cubic spline interpolation on a 2D grid as a function of $\epsilon$ and $z$. In model 2, it is calculated analytically as $n_{\text{EBL}}(\epsilon,z)=N_{0}(1+z)^{3.1}\epsilon^{-2.5}$ for $z\leq 1.4$ and $N_{0}(1+1.4)^{3.1}\epsilon^{-2.5}$ for larger $z$. In model 3, it is calculated as $n_{\text{EBL}}(\epsilon,z)=S(z)n_{0}(\epsilon/(1+z))$, with both factors found via cubic spline interpolation on 1D grids as functions of $z$ and $\epsilon/(1+z)$ respectively.

4.2 Photodisintegration cross sections

Measurements of photodisintegration cross sections are not available for all nuclides, and when there are, sometimes only the total and/or the one-neutron ejection cross sections have been measured. Cross sections for exclusive channels in which charged fragments are ejected are hard to measure, because such ejectiles tend to undergo multiple scattering in the target. Various phenomenological models have been used in UHECR propagation studies to treat these processes, with sizeable differences among them [40].

The following photodisintegration models are available in SimProp, and can be chosen via the -M command-line option:

[TABLE]

In photodisintegration models -M 0, -M 1 and -M 2, the processes considered are single nucleon ejection, double nucleon ejection, and multiple nucleon ejection; the type of each nucleon ejected is sampled at random (i.e., if a nucleus with 26 protons and 30 neutrons loses a nucleon, it is assumed to be a proton with probability $26/56$ and a neutron with probability $30/56$); this simplifying assumption is only approximately realistic (as in reality interaction channels yielding stable nuclei are more likely) and may result in a slight overestimate of the number of beta decays (and resulting neutrinos, which in any event are subdominant with respect to those from pion decay). In photodisintegration models -M 3 and -M 4, the processes implemented are single nucleon ejection (with the nucleon type chosen at random) and alpha-particle ejection.

The model used by default is that by Puget, Stecker and Bredekamp (1976) [22] as refined by Stecker and Salamon (1999) [23]. It does not distinguish between protons and neutrons, treating only beta-decay stable isobars for each mass number $2\leq A\leq 4$ and $9\leq A\leq 56$ (51 nuclides in total). The processes it models are single- and double-nucleon ejection for $\epsilon^{\prime}<\epsilon_{1}=30~{}\mathrm{MeV}$ and multiple-nucleon ejection with tabulated branching ratios for $\epsilon_{1}<\epsilon^{\prime}<\epsilon_{\max}=150~{}\mathrm{MeV}$. All other processes, e.g. those in which deuterons or alpha particles are ejected, are neglected in this model. An exception is beryllium-9, for which the only process modelled is fragmentation into one nucleon and two alpha particles.

If -M is used with a non-zero value, the first row of the standard input must contain the number of nuclear species in the model (excluding free nucleons), the energy $\epsilon_{1}$ and the energy $\epsilon_{\max}$; in the PSB model (-M 0) these are 51, $30~{}\mathrm{MeV}$ and $150~{}\mathrm{MeV}$ respectively. Then the standard input should include one line for each nuclide, giving cross section parameters as follows:

[TABLE]

All energies must be given in MeV and all cross sections must be given in millibarns. The parameter values used with the default model (-M 0) are listed in table 1. The corresponding cross sections are:

Models 0 and 1

[TABLE]

where

[TABLE]

branching ratios $b_{i}$ from table 2.

Model 2

[TABLE]

Model 3

[TABLE]

Model 4

[TABLE]

For $A\leq 4$, PSB cross sections are used, regardless of the photodisintegration model chosen; in the models where the channels implemented are $\sigma_{N}$ and $\sigma_{\alpha}$, for these nuclei we use $\sigma_{N}(\epsilon^{\prime})=\sum_{n}n\sigma_{n}^{\text{PSB}}(\epsilon^{\prime})$ and $\sigma_{\alpha}(\epsilon^{\prime})=0$.

The SimProp v2r4 package contains the files xsect_BreitWigner_TALYS-1.0.txt and xsect_BreitWigner_TALYS-1.6.txt intended for use with -M 2, xsect_BreitWigner2_TALYS-1.6.txt intended for use with -M 3, and xsect_Gauss2_TALYS-restored.txt intended for use with -M 4. The parameter values they contain were found from fits to the

results of TALYS [41], a program that can simulate nuclear reactions for a variety of projectile types and a wide range of projectile energies, computing cross sections for all exclusive channels, $\sigma_{n_{\mathrm{n}}n_{\mathrm{p}}n_{\rm d}n_{\rm t}n_{\rm h}n_{\alpha}}$ being the cross section for the channel in which $n_{\mathrm{n}}$ neutrons, $n_{\mathrm{p}}$ protons, $n_{\rm d}$ deuterons, $n_{\rm t}$ tritium nuclei, $n_{\rm h}$ helium-3 nuclei and $n_{\alpha}$ helium-4 nuclei are ejected.

In xsect_BreitWigner_TALYS-1.0.txt and xsect_BreitWigner_TALYS-1.6.txt, the one-nucleon ejection cross section $\sigma_{1}$ for $\epsilon_{1}^{\min}<\epsilon^{\prime}<\epsilon_{1}=30~{}\mathrm{MeV}$ was fitted to the sum of TALYS one-neutron and one-proton ejection cross sections $\sigma_{100000}+\sigma_{010000}$; the two-nucleon ejection cross section $\sigma_{2}$ for $\epsilon_{2}^{\min}<\epsilon^{\prime}<\epsilon_{1}=30~{}\mathrm{MeV}$ was fitted to the sum of TALYS two-neutron, one-neutron-one-proton, two-proton, and one-deuteron ejection cross sections $\sigma_{200000}+\sigma_{110000}+\sigma_{020000}+\sigma_{001000}$; and the total cross section for $\epsilon_{1}<\epsilon^{\prime}<\epsilon_{\max}=150~{}\mathrm{MeV}$ was fitted to the TALYS total cross section, while using the same branching ratios as in the PSB model. The default settings of TALYS-1.0 and TALYS-1.6 were used respectively.

In xsect_BreitWigner2_TALYS-1.6.txt and xsect_Gauss2_TALYS-restored.txt,the cross sections for nucleon and alpha-particle ejection $\sigma_{N}$ and$\sigma_{\alpha}$ were fitted to

[TABLE]

The rationale for not considering multiple nucleon ejection in these models is that, while the statistical uncertainties associated with the available $\ln A$ measurements are in principle small enough to distinguish protons from helium-4, they are too large to distinguish consecutive intermediate nuclei, e.g. carbon-12 from carbon-13; a detailed discussion about this is found in the SimProp v2r3 documentation [6] and in ref. [2].

As discussed in ref. [2], in the cases relevant to UHECR propagation, the results from released versions of TALYS used with their default settings have been found to be in much worse agreement with the measured total photodisintegration cross section data than those from the preliminary version used in ref. [42]. For this reason, in the file xsect_Gauss2_TALYS-restored.txt the parameters were fitted to results from TALYS-1.6 with settings restored to those used in the preliminary version, as described in ref. [2].

Also, all versions of TALYS largely overestimate the cross sections for channels in which alpha particles are emitted for the few nuclides for which measured data for these channels are available; these channels are neglected altogether in the PSB model. (The TALYS-1.6 settings that we restored to those used in the preliminary version only affect the total cross sections but not the branching ratios.)

5 Flux of secondary neutrinos and gamma rays

The main theoretical interest in secondary neutrinos and gamma rays is that this radiation carries information about sources of UHECRs not otherwise available. For instance, all protons reaching the Earth with $E>1~{}\mathrm{EeV}$ must originate from sources at $z<1$, no matter how much energy they started with (figure 2), and this limit is even stricter for other nuclei (especially light ones). Since the flux of extragalactic nuclei with $E\lesssim 1~{}\mathrm{EeV}$ may be suppressed by diffusion in intergalactic magnetic fields and/or contaminated by the high-energy tail of galactic cosmic rays, this means that all information about sources at $z>1$ is lost when looking at charged particles alone.

On the other hand, neutrinos produced in the decay of (photoproduced) pions can reach the Earth unaffected by their propagation (except for the redshift energy loss and flavour oscillations) even if produced up to $z\sim 10$. Neutrino energies range from a few $\mathrm{PeV}$ (in the case of pion production on EBL photons) to a few $\mathrm{EeV}$ (with CMB photons) [25].

Photons from pion decay and electrons/positrons from pair production and pion decay are produced with similar energies and undergo a cascade of pair production ($\gamma_{\text{HE}}+\gamma_{\text{backgr}}\to\mathrm{e}^{+}+\mathrm{e}^{-}$) and inverse Compton scattering ($\mathrm{e}^{\pm}+\gamma_{\text{backgr}}\to\mathrm{e}^{\pm}+\gamma_{\text{HE}}$) with very short interaction lengths, eventually resulting in cascades of photons with up to a few hundred GeV each. The shape of the cascade spectrum at Earth does not depend on the energy of the initiating partilce, provided it is sufficiently large $E_{\text{prod}}\gtrsim 100$ TeV [36]. Therefore, models of UHECR sources which only differ in the source emissivities at large $z$ and would be hard to distinguish from the observation of nuclei alone can still predict very different fluxes of neutrinos and gamma rays [5, 43, 10]. In principle neutrinos carry more information than the gamma rays, but they are harder to detect, therefore constraints from gamma-ray observations can be more stringent.

Fluxes of secondary gamma rays and neutrinos are also sensitive to the mass composition of UHECRs. For a given energy, heavy nuclei have a lower Lorentz factor and are less likely to undergo pion production with CMB photons. Therefore, the detection of $\mathrm{EeV}$ neutrinos or gamma rays would be evidence for a sizeable fraction of protons among the highest-energy cosmic rays.

In addition to this, whereas charged particles can be deflected by intergalactic and galactic magnetic fields, neutral particles travel in straight lines so that their arrival direction matches the angular position of the source.

5.1 Example scenarios

As an example of the dependence of gamma-ray fluxes on the assumptions about UHECR sources, in figure 3 we show the expected gamma-ray flux from UHECR propagation computed in the three pure-proton scenarios and the three mixed-composition scenarios of ref. [25]. The EBL model used is -L 7 (Gilmore et al. 2012 fiducial [29]). The source emissivities (table 3) were normalized to the latest Telescope Array surface detector data [44] in the proton-only scenario and the latest Pierre Auger Observatory combined spectrum data [45] above $10^{18}~{}\mathrm{eV}$ in the mixed-composition scenario, while all other parameters were kept fixed to the ref. [25] values.

The gamma-ray spectrum produced by cascading secondary gamma rays or pairs at each given redshift was computed using the analytic approximation of ref. [36] (their eqs. (9, 10), with $\mathcal{E}_{X}$ and $\mathcal{E}_{\gamma}$ interpolated as a function of redshift from the values shown in their figures 3 and 4).

Fermi-LAT measurements [46] of the isotropic diffuse gamma-ray background (IGRB) and the total extragalactic gamma-ray background (EGB) are also shown for comparison, the difference being the emissions resolved into point sources, which are excluded from the former and included in the latter. Note however that a sizeable fraction of the IGRB may originate in point sources that Fermi-LAT was unable to resolve, making the truly diffuse emission lower [9], but on the other hand if intergalactic magnetic fields are weak enough gamma-rays from UHCER propagation may come from the direction of the original source [47] and be excluded from the IGRB.

The model implemented here, based on the analytic approximation of [36], is very crude, neglecting the fact that the EBL is not monochromatic, and more realistic simulations of the cascade development would give more accurate results. The examples discussed here show how scenarios which would be practically indistinguishable from UHECR observations alone can result in gamma-ray fluxes differing by more than one order of magnitude, allowing observations to discriminate them. Results by other authors with slightly different assumptions [9, 48, 49] are qualitatively similar, but sometimes more stringent because they added to the predictions (or subtracted from the measurements) expected gamma-ray fluxes from other processes. A more extensive study of this topic will be the subject of a future dedicated publication.

Acknowledgments

We would like to thank CRPropa developers Rafael Alves Batista, Arjen van Vliet and David Walz, for helping us determine the origin of the differences between results from early versions of SimProp and CRPropa, and TALYS developers Arjan Koning and Stéphane Goriely, for their help in understanding the differences between the various versions of TALYS. We would also like to thank Pasquale Blasi, Carmelo Evoli and Oleg Kalashev for fruitful discussions about the electromagnetic cascades initiated by high-energy electrons and photons in intergalactic photon backgrounds.

The work by DB has received funding from the European Union’s Horizon 2020 research and innovation programme (Grant No. 646623). The work by AdM is supported by the IISN project 4.4502.16.

Appendix A Input parameters and output files

SimProp v2r4 recognizes the following command-line options:

The output is written in a ROOT file, whose name encodes the values of the parameters used (e.g. SimProp-v2r4_N100_A56_e17.0_E21.0_g1.00_z0.00_Z1.00_r0.00_L1_M0_n1.00_a1.00_D1_S1_p0_o0_s65539.root when the default settings are used). The file also contains the input parameter values in a tree called opt.

If the parameter -o is set to 0, the file contains the following trees:

**Tree **nuc

This tree has an entry for each branch. The name of the tree is due to historical reasons; particles other than nuclei are also included now.

**Tree **ev

This tree has an entry for each event.

If the parameter -o is set to 1, the file contains the following tree:

**Tree **summary

This tree has an entry for each event.

If the parameter -o is set to 2, the file contains all the trees described above.

Appendix B Distance measures

In an expanding universe, there are several possible definitions of distance which for sizeable $z$ are not equivalent. The following definitions are valid for flat space ($\Omega=\Omega_{\text{m}}+\Omega_{\Lambda}=1$).

The comoving distance is the proper distance between the positions of two objects measured at a fixed time, divided by the scale factor $a(t)=R(t)/R_{0}=(1+z)^{-1}$ at that time. The comoving distance of two objects moving with the Hubble flow does not vary with time. The comoving distance of an object whose light reaches us today after leaving the object at redshift $z$ is given by

[TABLE]

Likewise, the comoving volume is given by the proper volume times $(1+z)^{3}$; the number density per unit comoving volume of a fixed number of objects moving with the Hubble flow does not vary with time.

The light travel distance is the cosmological time elapsed since the light leaves an object until it reaches us. It is given by

[TABLE]

In the case of a distribution of closely spaced identical sources, we can define the source emissivity $\mathcal{L}$ as the total energy injected per unit comoving volume per unit time, i.e. $\mathcal{L}=n_{\text{s}}L$, where $n_{\text{s}}$ is the number density of sources per unit comoving volume, and the luminosity $L$ of each source is the total energy emitted by the source per unit time,

[TABLE]

$Q$ being the injection spectrum (number of particles emitted per unit energy per unit time) of each source. Likewise, we can define $\mathcal{Q}(E_{\text{inj}})=n_{\text{s}}Q(E_{\text{inj}})$.

In the cases where the injection spectrum and density of sources depends on energy and possibly time (or equivalently redshift) but not on position, e.g. $\mathcal{Q}=\mathcal{Q}(E_{\text{inj}},t)$, it can be shown that the expected fluxes at Earth are given by

[TABLE]

where $T_{ij}(E_{\text{Earth}}|E_{\text{inj}},t)$ is the average number of particles of type $i$ with energy $E_{\text{Earth}}$ at the present time $t_{0}$ from each particle of type $j$ injected with energy $E_{\text{inj}}$ at time $t$. Therefore, in the case $\mathcal{Q}(E_{\text{inj}},t)=\mathcal{Q}_{0}(E_{\text{inj}})S(z)$, in order to correctly analyse simulations of UHECR propagation in which the source positions are sampled from a uniform distribution in $z$, each event must be weighed by a factor proportional to $S(z)\absolutevalue{\differential{t}/\differential{z}}$.

Appendix C Performances

In table 8, we list the CPU times used by SimProp v2r4 simulations in a few example configurations on machines with AMD Opteron 6238 CPUs with 6 cores and 2.6 GHz clock rates. The times mainly depend on the choice of settings about the production of secondary neutrinos, electrons and photons: producing none (-S 0 -p 0) is usually two orders of magnitude faster than producing all, binning electrons and photons by redshift (-S 2 -p 2). All other things being equal, the computation time is very roughly proportional to the injection mass number and redshift, whereas the injection energy does not have a major impact.