SIMPLIFIED MODELS FOR PHOTOHADRONIC INTERACTIONS IN COSMIC ACCELERATORS

Photohadronic interactions in high-energy astrophysical accelerators are the key ingredient of hadronic source models besides proton–proton interactions. The smoking gun signature of these interactions may be the neutrino production from charged pions, which could be detected in neutrino telescopes (ANTARES Collaboration 1999; Karle et al. 2003; Tzamarias 2003; NEMO Collaboration 2005), such as IceCube; see, e.g., Rachen & Mészáros (1998) for the general theory of the astrophysical sources. For instance, in gamma-ray bursts (GRBs), photohadronic interactions are expected to lead to a significant flux of neutrinos in the fireball scenario (Waxman & Bahcall 1997). On the other hand, in active galactic nucleus (AGN) models, the neutral pions produced in these interactions may describe the second hump in the observed photon spectrum, depending on the dominance of synchrotron or inverse Compton cooling of the electrons. The protons in these models are typically assumed to be accelerated in the relativistic outflow together with electrons and positrons by Fermi shock acceleration. The target photon field is typically assumed to be the synchrotron photon field of the co-accelerated electrons and positrons. Also thermal photons from broad-line regions, the accretion disk, and the cosmic microwave background (CMB) may serve as the target photon field. The latter case, relevant for the cosmogenic neutrino flux, is not considered in detail.

The basic ideas of complete hadronic models for AGNs have already been described in previous works (Mannheim 1993; Mücke & Protheroe 2001; Aharonian 2002), and leptonic models have been discussed by Maraschi et al. (1992) and Dermer & Schlickeiser (1993). In the first place, these models have been used as static models to describe steady-state spectral energy distributions, but with today's generation of gamma-ray telescopes, a detailed analysis of the dynamics of very high energetic sources is possible, and time-dependent modeling is inevitable. For the case of leptonic models, see, e.g., Böttcher & Chiang (2002) and Rüger et al. (2010). A necessary prerequisite for a time-dependent hadronic modeling is the efficient computation of the photohadronic interactions. The online calculation via Monte Carlo simulations is not feasible, therefore, a parametric description is the most viable method; see, e.g., Kelner & Aharonian (2008).

The prediction of neutrino fluxes in many source models relies on the photohadronic neutrino production. In this case, astrophysical neutrinos are normally assumed to originate from pion decays, with a flavor ratio at the source of (f_e,  f_μ,  f_τ) ≃ (1/3,  2/3,  0) arising from the decays of both primary pions and secondary muons. However, it was pointed out in Rachen & Mészáros (1998) and Kashti & Waxman (2005) that such sources may become opaque to muons at higher energies, in which case the flavor ratio at the source changes to (f_e,  f_μ,  f_τ) ≃ (0,  1,  0). Therefore, one can expect a smooth transition from one type of source to the other as a function of the neutrino energy (Kachelrieß & Tomàs 2006; Lipari et al. 2007; Kachelrieß et al. 2008), depending on the cooling processes of the intermediate muons, pions, and kaons. Recently, the use of flavor information has been especially proposed to extract some information on the particle physics properties of the neutrinos and the properties of the source; see Pakvasa (2008) for a review. For instance, if the neutrino telescope has some flavor identification capability, this property can be used to extract information on the decay (Beacom et al. 2003; Lipari et al. 2007; Majumdar 2007; Maltoni & Winter 2008; Bhattacharya et al. 2009) and oscillation (Farzan & Smirnov 2002; Beacom et al. 2004; Serpico & Kachelrieß 2005; Serpico 2006; Bhattacharjee & Gupta 2005; Winter 2006; Majumdar & Ghosal 2007; Meloni & Ohlsson 2007; Blum et al. 2007; Rodejohann 2007; Xing 2006; Pakvasa et al. 2008; Hwang & Siyeon 2008; Choubey et al. 2008; Esmaili & Farzan 2009) parameters, in a way which might be synergistic to terrestrial measurements. Of course, the flavor ratios may also be used for source identification; see, e.g., Xing & Zhou (2006) and Choubey & Rodejohann (2009). Except for flavor identification, the differentiation between neutrinos and antineutrinos could be useful for the discrimination between pγ and pp induced neutrino fluxes, or for the test of neutrino properties (see, e.g., Maltoni & Winter 2008). A useful observable may be the Glashow resonance process, $\bar{\nu }_e + e^-\rightarrow W^-\rightarrow {\rm anything}$, at around 6.3 PeV (Learned & Pakvasa 1995; Anchordoqui et al. 2005; Bhattacharjee & Gupta 2005) to distinguish between neutrinos and antineutrinos in the detector. Within the photohadronic interactions, the π⁺-to-π⁻ ratio determines the ratio between electron neutrinos and antineutrinos at the source. Therefore, a useful source model for these applications should include sufficiently accurate flavor ratio predictions, including the possibility to include the cooling of the intermediate particles, as well as π⁺-to-π⁻ ratio predictions. In addition, the computation of the neutrino fluxes should be efficient enough to allow for reasonable parameter studies or to be used as a fit model.

Photohadronic interactions in astrophysical sources are typically either described by the refined Monte Carlo simulation of the SOPHIA software (Mücke et al. 2000a), which is partially based on Rachen (1996), or are in very simplified approaches computed with the Δ-resonance approximation

$$\begin{equation} p + \gamma \rightarrow \Delta ^+ \rightarrow \left\lbrace \begin{array}{@{}ll@{}} n + \pi ^+ & 1/3 \, \, {\rm of\ all\ cases} \\ p + \pi ^0 & 2/3 \, \, {\rm of\ all\ cases}. \\ \end{array} \right. \end{equation} \tag{ 1 }$$

The SOPHIA software probably provides the best state-of-the-art implementation of the photo-meson production, including not only the Δ-resonance, but also higher resonances, multi-pion production, and direct (t-channel) production of pions. Kaon production is included as well by the simulation of the corresponding quantum chromodynamics (QCD) processes (fragmentation of color strings). The treatment in Equation (1), on the other hand, has been considered sufficient for many purposes, such as estimates for the neutrino fluxes. However, both of these approaches have disadvantages: the statistical Monte Carlo approach in SOPHIA is too slow for efficient use in every step of modern time-dependent source simulations of AGNs and GRBs. The treatment in Equation (1), on the other hand, does not allow for predictions of the neutrino–antineutrino ratio and the shape of the secondary particle spectra because higher resonances and other processes contribute significantly to these. One possibility for obtaining a more accurate model is to use different interaction types with different (energy-dependent) cross sections and inelasticities (fractional energy loss of the initial nucleon). For example, one may define an interaction type for the resonances, and an interaction type for multi-pion production (see, e.g., Reynoso & Romero 2009; Mannheim & Biermann 1989, and others). Typically, such approaches do not distinguish π⁺ from π⁻ production. An interesting alternative has been proposed in Kelner & Aharonian (2008), which approximates the SOPHIA treatment analytically. It provides a simple and efficient way to compute the electron, photon, and neutrino spectra by integrating out the intermediate particles. Naturally, the cooling of the intermediate particles cannot be included in such an approach.

Since we are also interested in the cooling of the intermediate particles (muons, pions, kaons), we follow a different direction. We propose a simplified model based on the very first physics principles, which means that the underlying interaction model can be easily adapted if new data are provided. In this approach, we explicitly include the intermediate particles. We illustrate the results for the neutrino spectra, but the extension to photons and electrons/positrons is straightforward.

More explicitly, the requirements for our interaction model are the following.

• 1.  
  The model should predict the π⁺, π⁻, and π⁰ fluxes separately, which is needed for the prediction of photon, neutrino, and antineutrino fluxes, and their ratios.
• 2.  
  The model should be fast enough for time-dependent calculations and for systematic parameter studies. This, of course, requires compromises. For example, compared to SOPHIA, our kinematics treatment will be much simpler.
• 3.  
  The particle physics properties should be transparent, easily adjustable, and extendable.
• 4.  
  The model should be applicable to arbitrary proton and photon input spectra.
• 5.  
  The secondaries (pions, muons, kaons) should not be integrated out because (a) their synchrotron emission may contribute to observations, (b) the muon (and pion/kaon) cooling affects the flavor ratios of neutrinos, and (c) pion cooling may be in charge of a second spectral break in the prompt GRB neutrino spectrum; see Waxman & Bahcall (1999).
• 6.  
  The cooling and escape timescales of the photohadronic interactions as well as the proton/neutron re-injection rates should be provided, since these are needed for time-dependent and steady-state models.
• 7.  
  The kaons leading to high-energy neutrinos should be roughly provided, since their different cooling properties may lead to changes in the neutrino flavor ratios (see, e.g., Kachelrieß & Tomàs 2006; Lipari et al. 2007; Kachelrieß et al. 2008). Therefore, we incorporate a simplified kaon production treatment to allow for the test of the impact of such effects. This also serves as a test case for how to include new processes.

For the underlying physics, we mostly follow similar principles as in SOPHIA.

This work is organized as follows. In Section 2, we review the basic principles of the particle interactions, in particular, the necessary information (response function) to compute the meson photoproduction for arbitrary photon and proton input spectra. In Section 3, we summarize the meson photoproduction as implemented in the SOPHIA software, but we simplify the kinematics treatment. In Section 4, we define simplified models based on that. As a key component, we define appropriate interaction types such that we can factorize the two-dimensional response function. In Section 5, we then summarize the weak decays and compare the different approaches with each other and with the output from the SOPHIA software. For the sake of completeness, we provide in Appendix B how the cooling and escape timescales and the neutron/proton re-injection can be computed in our simplified models.

For the notation, we follow Lipari et al. (2007), where we first focus on weak decays, such as π⁺ → μ⁺ + ν_μ, and then extend the discussion to photohadronic interactions. Let us first consider a single decay process. The production rate Q_b(E_b) (per energy interval) of daughter particles of species b and energy E_b from the decay of the parent particle a can be written as

$$\begin{equation} Q_b(E_b) = \sum _{\mathrm{IT}} \int dE_a \, N_a(E_a) \, \Gamma ^{\mathrm{IT}}_{a \rightarrow b}(E_a) \, \frac{d n_{a \rightarrow b}^{\mathrm{IT}}}{dE_b} (E_a,E_b) \,. \end{equation} \tag{ 2 }$$

Here, N_a(E_a) is the differential spectrum of parent particles¹ (particles per energy interval), and Γ^IT_a→b is the interaction or decay rate (probability per unit time and particle) for the process a → b as a function of energy E_a (which is assumed to be zero below the threshold). Since in pion photoproduction many interaction types contribute, we split the production probability in interaction types "IT."

The function dn^IT_a→b/dE_b(E_a, E_b) describes the distribution (as a function of parent and daughter energy) of daughter particles of type b per final state energy interval dE_b. This function can be non-trivial. It contains the kinematics of the decay process, i.e., the energy distribution of the discussed daughter particle. Other species, which we are not interested in, are typically integrated out. If more than one daughter particle of the same species b is produced, or less than one (in average) because of other branchings, it must also give the number of daughter particles per event as a function of energy, which is often called "multiplicity."

Note that the decay can be calculated in different frames, such as the parent rest frame (PRF), in the center of mass frame of the parents (CMF), or in the shock rest frame (SRF), typically used to describe shock accelerated particles (such as from Fermi shock acceleration) in astrophysical environments. However, the cross sections, entering the interaction rate Γ, are often given in a particular frame, which has to be properly included. In addition, note that Q and N are typically given per volume, but here this choice is arbitrary since it enters on both sides of Equation (2).

Sometimes it is useful to consider the interaction or decay chain a → b → c without being interested in b, for instance, for the decay chain π → μ → ν. In this case, one can integrate out b. This approach has, for instance, been used in Kelner & Aharonian (2008) for obtaining the neutrino spectra. Note, however, that the parent particles a or b may lose energy before they interact or decay, such as by synchrotron radiation. We do not explicitly treat such energy losses in this study, but we provide a framework to include them, such as in a steady-state or time-dependent model for each particle species.

In meson photoproduction, a similar mechanism can be used. The interaction rate (Equation (2)) can be interpreted in terms of the incident protons (or neutrons) because of the much higher energies in the SRF, i.e., the species a is identified with the proton or neutron, which we further abbreviate as p or p' (for proton or neutron). In this case, the interaction rate depends on the interaction partner, the photon, as

$$\begin{eqnarray} \Gamma ^{\mathrm{IT}}_{p \, \gamma \rightarrow p^{\prime } \, b}(E_p) &=& \int d \varepsilon \int _{-1}^{+1} \frac{d \cos \theta _{p \gamma }}{2} \, (1- \cos \theta _{p \gamma }) \nonumber\\ &&\times n_\gamma (\varepsilon,\cos \theta _{p \gamma }) \, \sigma ^{\mathrm{IT}}(\epsilon _r) \,. \end{eqnarray} \tag{ 3 }$$

Here, n_γ(ε, cos θ_pγ) is the photon density as a function of photon energy ε and the angle between the photon and proton momenta θ_pγ, σ^IT(_r) is the photon production cross section, and _r = E_pε/m_p(1 − cos_pγ) is the photon energy in the nucleon/PRF in the limit β_p ≈ 1. The interaction itself, and therefore E_p and ε, is described in the SRF. The daughter particles b are typically π⁺, π⁻, π⁰, or kaons. If intermediate resonances are produced, we integrate them out so that only pions (kaons) and protons or neutrons remain as the final states.

In the following, we will assume isotropy n_γ(ε, cos θ_pγ) ≃ n_γ(ε) of the photon distribution. This limits this specific model to scenarios where we have seed photons produced in the SRF (i.e., synchrotron emission). The other interesting case, where thermal photons are coming from outside the shock, may be easily implemented as well when the shock speed is high enough to assume delta-peaked angular distributions. Arbitrary angular distributions would require additional integrations. Assuming that the photon distribution is isotropic, the integral over cos θ_pγ can be replaced by one over _r, we have

$$\begin{equation} \Gamma ^{\mathrm{IT}}_{p \, \gamma \rightarrow p^{\prime } \, b}(E_p) = \frac{1}{2} \frac{m_p^2}{E_p^2} \int _{\frac{\epsilon _{\mathrm{th}} m_p}{2 E_p}}^{\infty } d \varepsilon \frac{n_\gamma (\varepsilon)}{\varepsilon ^2} \int _{\epsilon _{\mathrm{th}}}^{2 E_p \varepsilon /m_p} d \epsilon _r \, \epsilon _r \, \sigma ^{\mathrm{IT}}(\epsilon _r) \,. \end{equation} \tag{ 4 }$$

Here, _th ≃ 150 MeV is the threshold below which the cross sections are zero. Note that _r corresponds to the available center-of-mass energy $\sqrt{s}$ of the interaction as

$$\begin{equation} s(\epsilon _r) = m_p^2 + 2 \, m_p \, \epsilon _r \,. \end{equation} \tag{ 5 }$$

In Equation (4), the integral over _r takes into account that the proton and photon may hit each other in different directions. In the most optimistic case, _r ≈ 2E_pε/m_p (head-on hit), whereas in the most pessimistic case, _r ≈ 0 (photon and proton travel in the same direction). It should be noted that the assumption of isotropy here is limiting the model to internally produced photon fields. This could be lifted for different angular distributions, but except for the case of unidirectional photons, this would require one more integration.

In general, the function dn^IT_a→b/dE_b(E_a, E_b) in Equation (2) is a non-trivial function of two variables E_a and E_b. For photo-meson production in the SRF, the energy of the target photons is typically much smaller than the incident nucleon energy, and β_a ≃ 1. In this case,²

$$\begin{equation} \frac{d n_{a \rightarrow b}^{\mathrm{IT}}}{dE_b} (E_a,E_b) \simeq \delta \big(E_b - \chi _{a \rightarrow b}^{\mathrm{IT}} \, E_a\big) \cdot M_b^{\mathrm{IT}} \,. \end{equation} \tag{ 6 }$$

The function χ^IT_a→b, which depends on the kinematics of the process, describes which (mean) fraction of the parent energy is deposited in the daughter species. The function M^IT_b describes how many daughter particles are produced at this energy in average. For our purposes, it will typically be a constant number which depends on interaction type and species b (if it changes at a certain threshold, one can define different interaction types below and above the threshold). For the relationship between χ^IT_a→b and the inelasticity K (fractional energy loss of the initial nucleon), see Appendix B. As we will discuss later, Equation (6) is an oversimplified for more realistic kinematics, such as in direct or multi-pion production, since χ^IT_a→b is, in general, a more complicated function of E_b/E_a and the initial state properties. In this case, the distribution is not sufficiently peaked around its mean, and the δ-distribution in Equation (6) has to be replaced by a broader distribution function describing the distribution of the daughter energies. Instead of choosing a broader distribution function at the expense of efficiency, in this case we will define different interaction types with different values of χ^IT_a→b, simulating the broad energy distribution after the integration over the input spectra.

As the next step, we insert Equations (6) and (4), valid for photoproduction of pions, into Equation (2), in order to obtain

$$\begin{equation} Q_b(E_b) = \int _{E_b}^{\infty } \frac{dE_p}{E_p} \, N_p(E_p) \, \int _{\frac{\epsilon _{\mathrm{th}} m_p}{2 E_p}}^{\infty } d\varepsilon \, n_\gamma (\varepsilon) R_b(x,y), \end{equation} \tag{ 7 }$$

with

$$\begin{equation} x \equiv \frac{E_b}{E_p} \,, \quad y\equiv \frac{E_p \, \varepsilon }{m_p} \,, \end{equation} \tag{ 8 }$$

and the "response function"

$$\begin{eqnarray} R_b(x,y) &\equiv &\sum _{\mathrm{IT}} R^{\mathrm{IT}}(x,y)\nonumber\\ & \equiv & \sum _{\mathrm{IT}} \frac{1}{2 y^2} \! \int _{\epsilon _{\mathrm{th}}}^{2 y} \! d \epsilon _r \, \epsilon _r \, \sigma ^{\mathrm{IT}}(\epsilon _r) M^{\mathrm{IT}}_b(\epsilon _r) \, \delta (x - \chi ^{\mathrm{IT}}(\epsilon _r)).\nonumber\\ \end{eqnarray} \tag{ 9 }$$

Here, 0 < x = E_b/E_p < 1 is the secondary energy as a fraction of the incident proton or neutron energy, _th ≲ 2y = 2E_pε/m_p ≲ 10⁴ GeV corresponds to the maximal available center-of-mass energy,³ and the χ^IT is the fraction of proton energy deposited in the secondary. Note that χ^IT(_r) and M_b(_r) are typically functions of _r, if they only depend on the center-of-mass energy of the interaction. We will define our interaction types such that M_b is independent of _r.

If the response function in Equation (9) is known, the secondary spectra can be calculated from Equation (7) for arbitrary injection and photon spectra. A similar approach was used in Kelner & Aharonian (2008) for the neutrinos and photons directly. However, we do not integrate out the intermediate particles, which in fact leads to a rather complicated function R(x, y), summed over various interaction types. In Figure 1 we show the cross section as a function of _r for these interaction types separately, where the baryon resonances have been summed up. Naively, one would just choose the dominating contributions in the respective energy ranges in order to obtain a good approximation for σ^IT(_r). However, the different contributions have different characteristics, such as different π⁺-to-π⁻ ratios in the final states, and therefore different neutrino–antineutrino ratios. In addition, the function δ(x − χ^IT(_r)) in Equation (9) maps the same region in _r in different regions of x = E_b/E_p, depending on the interaction type. This means that while a particular cross section may dominate for certain energies _r, for instance, the pion energies where the interaction products are found can be very different, leading to distinctive features. Therefore, for our purposes (such as flavor ratio and neutrino–antineutrino predictions), it is not a sufficient approximation to just choose the dominating cross section.

Zoom In Zoom Out Reset image size

From the particle physics point of view, R_b(x, y) = R_b(E_b, E_p, ε) is very similar to a detector response function in a fixed target experiment. It describes the reconstructed particle energy spectrum as a function of the properties of the incident proton beam. As the major difference, the "detector material" is kept as a variable function of energy, leading to the second integral over the photon density.

The description of photohadronic interactions is based on Rachen (1996) and Mücke et al. (2000a), which means that the fundamental physics are similar to SOPHIA. However, our kinematics will be somewhat simplified. The purpose of this section is to show the key features of the interaction model. In addition, it is the first step toward an analytical description for the response function in Equation (9), which should be as simple as possible for our purposes, and as accurate as necessary. Note that we do not distinguish between protons and neutrons for the cross sections (there are some small differences in the resonances and the multi-pion production).

3.1. Summary of Processes

In summary, we consider the following processes.

${ \bDelta}$-resonance region. The dominant resonance process is the Δ(1232)-resonance (at _r = 340 MeV; cf. Equation (5) for the _r–s conversion):

$$\begin{equation} p + \gamma \stackrel{\Delta (1232)}{\rightarrow } p^{\prime } + \pi \,. \end{equation} \tag{ 10 }$$

This process produces neutral (for p' = p) or charged (for p' ≠ p) pions (e.g., for protons in the initial state), see Equation (1).

Higher resonances. The most important higher resonance contribution is the decay chain

$$\begin{equation} \gamma + p \stackrel{\Delta, N}{\rightarrow } \Delta ^{\prime } + \pi \,, \quad \Delta ^{\prime } \rightarrow p^{\prime } + \pi ^{\prime } \end{equation} \tag{ 11 }$$

via Δ- and N-resonances. Other contributions, we consider, come from the decay chain

$$\begin{equation} \gamma + p \stackrel{\Delta, N}{\rightarrow } \rho + p^{\prime } \,, \quad \rho \rightarrow \pi + \pi ^{\prime } \,. \end{equation} \tag{ 12 }$$

Direct production. The same interactions as in Equation (10) or Equation (11) (with the same initial and final states) can also take place in the t-channel, meaning that the initial γ and nucleon exchange a pion instead of creating a (virtual) baryon resonance in the s-channel, which again decays. This mechanism is also called "direct production" because the properties of the pion are already determined at the nucleon vertex. For instance, the process p + γ → p + π⁰ only takes place in the s-channel, whereas in the process p + γ → n + π⁺ the t-channel is possible as well, because the photon can only couple to charged pions. The s-channel reactions typically have a pronounced peak in the _r-distribution, whereas they are almost structureless in the momentum transfer distribution. On the other hand, t-channel reactions do not have the strong peak in the center-of-mass energy, but have a strong correlation between the initial and final state momentum distributions, implying that the pions are found at lower energies. On a logarithmic scale, however, the momentum distribution of the pions is broad.

Multi-pion production. At high energies, the dominant channel is statistical multi-pion production leading to two or more pions. In SOPHIA (Mücke et al. 2000a) the process is described by the QCD fragmentation of color strings.

The effect of kaon decays is usually small. However, kaon decays may have interesting consequences for the neutrino flavor ratios at very high energies, in particular, if strong magnetic fields are present (the kaons are less sensitive to synchrotron losses because of the larger rest mass; Kachelrieß & Tomàs 2006; Kachelrieß et al. 2008). Therefore, we consider the leading mode: K⁺ production (for protons in the initial state) with the decay channel leading to highest energy neutrinos.⁴ Note that at even higher energies, other processes, such as charmed meson production, may contribute as well.

In Figure 1 we show the total pγ photo-meson cross section as a function of _r, analogously to Mücke et al. (2000a) (1μbarn = 10⁻³⁰ cm²; data from Particle Data Group 2008). The contributions of baryon resonances (red, dashed), the direct channel (green, dotted), and multi-pion production (brown) are shown separately.

In order to fully describe Equation (9) for each interaction type, we need the kinematics, entering χ^IT, the multiplicities M^IT_b, and the cross section σ^IT.

3.2. Kinematics and Secondary Multiplicities

The kinematics for the resonances and direct production can be effectively described in the two-body picture. We follow the calculation of Berezinsky & Gazizov (1993) for the reaction p + γ → p' + π with E_p ≫ E_γ (β_p ≈ 1) in the SRF. The Lorentz factor between CMF and SRF is

$$\begin{equation} \gamma _{\rm cm}=\frac{E_p+E_\gamma }{\sqrt{s}}\simeq \frac{E_p}{\sqrt{s}}, \end{equation} \tag{ 13 }$$

with $\sqrt{s}$ being the total CMF energy from Equation (5). The energy of the pion in the SRF can then be written as

$$\begin{equation} E_\pi =\gamma _{\rm cm} E_\pi ^{\rm cm} \big(1+ \beta _\pi ^{\rm cm}\cos {\theta _\pi ^{\rm cm}}\big) = E_p \, \frac{E_\pi ^{\rm cm}}{\sqrt{s}}\big(1+\beta _\pi ^{\rm cm}\cos {\theta _\pi ^{\rm cm}}\big), \end{equation} \tag{ 14 }$$

with E^cm_π, p^cm_π, and θ^cm_π, which are the pion energy, momentum, and angle of emission in the CMF, respectively. Note that backward scattering of the pion means that cos θ^cm_π = −1. From Equation (14), we can read off the fraction of energy of the initial p going into the pion E_π/E_p (which is equal to the inelasticity of p for p' = p). Since p + γ → p' + π is a two-body reaction, the energies of p' and pion are given by (Particle Data Group 2008)

$$\begin{equation} E_{\pi }^{\rm cm}=\frac{s-m_p^2+m_\pi ^2}{2\sqrt{s}}, \end{equation} \tag{ 15 }$$

with s(_r) given by Equation (5). Using Equation (15) in Equation (14), we can calculate χ in Equation (9) as

$$\begin{equation} \chi (\epsilon _r)= \frac{s(\epsilon _r)-m_p^2+m_\pi ^2}{2 s(\epsilon _r)} \big(1+\beta _\pi ^{\mathrm{cm}} \cos {\theta _\pi }\big) \,. \end{equation} \tag{ 16 }$$

Indeed, χ(_r) is a function of _r if cos θ_π is a constant. As described in Rachen (1996), the average <cos θ_π> ≃ 0 for the resonances (such as for isotropic emission in the CMF), whereas the direct production is backward peaked to a first approximation for sufficiently high energies <cos θ_π> → −1. This approximation is only very crude. Therefore, we calculate the mean <cos θ_π> for direct production by averaging the probability distribution of the Mandelstam variable t, as we discuss in Appendix A. Note that even for the resonances these scattering angle averages are only approximations; a more refined kinematics treatment, such as in SOPHIA, will lead to a smearing around these mean values. In addition, note that we also use Equation (16) for the kaon production (with the rough estimate <cos θ_π> ≃ 0), where the replacements m_p → m_Λ and m_π → m_K have to be made. Similarly, one has for the first and second pion for the interaction in Equation (11) (Rachen 1996)

$$\begin{equation} \chi _a(\epsilon _r)=\frac{s(\epsilon _r)-m_\Delta ^2+m_\pi ^2}{2 s(\epsilon _r)} \big(1+\beta _\pi ^{{\rm cm}} \cos {\theta _\pi }\big) \,, \end{equation} \tag{ 17 }$$

$$\begin{equation} \chi _b(\epsilon _r)=\frac{s(\epsilon _r)-m_\pi ^2+m_\Delta ^2}{2 s(\epsilon _r)} \frac{m_\Delta ^2-m_p^2+m_\pi ^2}{2 m_\Delta ^2} \big(1+\beta _\Delta ^{{\rm cm}} \cos {\theta _\Delta }\big) \,, \end{equation} \tag{ 18 }$$

where Δ = Δ(1232), and for the interaction Equation (12)

$$\begin{equation} \chi (\epsilon _r)=\frac{1}{2}\frac{s(\epsilon _r)-m_p^2+m_\rho ^2}{2 s(\epsilon _r)} \big(1+\beta _\rho ^{{\rm cm}} \cos {\theta _\rho }\big) = \frac{K}{2}, \end{equation} \tag{ 19 }$$

with m_ρ ≃ 775 MeV. Note that in order to evaluate the δ-function in Equation (9) (if one integrates over _r), one also needs the derivative χ'(_r), which can easily be computed from the functions above.

The average multiplicities for different types of resonances and direct production channels can be found in Table 1. The multiplicities of the pions M_π add up to the number of pions produced (one or two), the multiplicities for the nucleons M'_p in the final state p' to one. They are needed in Equations (6) and (B1) (cooling and escape timescale in Appendix B). In the last columns, the cosine of the approximate scattering angle in the center of mass frame is given, together with the equation for kinematics. The interaction types IT are labeled by the type of resonance (Δ or N) or T for the direct production (t-channel). Note that different resonances will contribute a similar pattern, such as through the interaction types Δ₁(1232), Δ₁(1700), Δ₁(1905), Δ₁(1950), etc. In addition note that the branching ratios are absorbed in the cross sections. Therefore, the total yield is always one (or two) pions per interaction in Table 1.

Table 1. Summary of the Average Multiplicities for Different Types of Considered Resonances and Direct Production Channels

IT	Initial p is Proton					Initial p is Neutron					Kinematics
	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	Mn	Mp	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	Mn	Mp	Equation	<cos θ>
Interactions $\gamma + p \stackrel{\mathrm{IT}}{\rightarrow } p^{\prime } + \pi$
Δ1	$\frac{2}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	(16)	0
N1	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	(16)	0
T1	⋅⋅⋅	1	⋅⋅⋅	1	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	1	⋅⋅⋅	1	(16)	−1
Interactions $\gamma + p \stackrel{\mathrm{IT}}{\rightarrow } \Delta ^{\prime } + \pi$, Δ' → p' + π'
Properties of first pion π and nucleon p'
Δ2a	$\frac{1}{15}$	$\frac{8}{15}$	$\frac{2}{5}$	$\frac{17}{45}$	$\frac{28}{45}$	$\frac{1}{15}$	$\frac{2}{5}$	$\frac{8}{15}$	$\frac{28}{45}$	$\frac{17}{45}$	(17)	0
N2a	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{2}$	$\frac{2}{9}$	$\frac{7}{9}$	$\frac{1}{3}$	$\frac{1}{2}$	$\frac{1}{6}$	$\frac{7}{9}$	$\frac{2}{9}$	(17)	0
T2a	⋅⋅⋅	$\frac{1}{4}$	$\frac{3}{4}$	$\frac{1}{6}$	$\frac{5}{6}$	⋅⋅⋅	$\frac{3}{4}$	$\frac{1}{4}$	$\frac{5}{6}$	$\frac{1}{6}$	(17)	−1
Properties of second pion π'
Δ2b	$\frac{2}{5}$	$\frac{19}{45}$	$\frac{8}{45}$			$\frac{2}{5}$	$\frac{8}{45}$	$\frac{19}{45}$			(18)	0
N2b	$\frac{1}{3}$	$\frac{11}{18}$	$\frac{1}{18}$			$\frac{1}{3}$	$\frac{1}{18}$	$\frac{11}{18}$			(18)	0
T2b	$\frac{1}{6}$	$\frac{3}{4}$	$\frac{1}{12}$			$\frac{1}{6}$	$\frac{1}{12}$	$\frac{3}{4}$			(18)	1
Interactions $\gamma + p \stackrel{\mathrm{IT}}{\rightarrow } \rho + p^{\prime }$, ρ → π + π' (π and π' summed over)
Δ3	$\frac{1}{3}$	1	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	1	$\frac{2}{3}$	$\frac{1}{3}$	(19)	0
N3	$\frac{2}{3}$	1	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	1	$\frac{1}{3}$	$\frac{2}{3}$	(19)	0

Notes. The multiplicities of the pions M_π add up to the number of pions produced (one or two), the multiplicities for the nucleons M'_p in the final state p' to one. In the last columns, the cosine of the approximate scattering angle in the center of mass frame is given, together with the equation for kinematics. The interaction types IT are labeled by the type of resonance (Δ or N) or "Direct" for the direct production (t-channel).

Download table as:  ASCIITypeset image

For the multi-pion channel, we use two different approaches. A very simple treatment can be performed similar to Atoyan & Dermer (2003), with some elements of Mücke et al. (2000a). Most of the energy lost by the proton (≃0.6 E_p or K ≃ 0.6) is split equally among three pions. The three types π⁺, π⁻, and π⁰ are therefore approximately produced in equal numbers. Our multi-pion channel, parameterized by the multi-pion cross section in Figure 1, however, is actually a combination of different processes and residual cross sections. For instance, diffractive scattering is a small contribution. Therefore, in order to reproduce the pion multiplicity time cross sections in Figures 9 and 10 of Mücke et al. (2000a) more accurately, we choose $M_{\pi ^0}=1$ ($M_{\pi ^0}=1$), $M_{\pi ^+}=1.2$ ($M_{\pi ^+}=0.85$), and $M_{\pi ^-} = 0.85$ ($M_{\pi ^-} =1.2$) for initial protons (neutrons). Close to the threshold, the decreasing phase space for pion production requires a modification of the threshold _r ⩾ 1 GeV for π⁻ (π⁺) from initial protons (neutrons). This corresponds to a vanishing multiplicity below the threshold, i.e., we assume that below _r ⩾ 1 GeV only two pions are produced. We assume for the fraction of the proton energy going into one pion produced in the multi-pion channel χ^Multi-π ≃ 0.2. In addition, we choose M_p = 0.69 and M_n = 0.31 for initial protons (or M_p = 0.31 and M_n = 0.69 for initial neutrons) in accordance with Figure 11 of Mücke et al. (2000a) for high energies. As an alternative, we use a more accurate but computationally more expensive approach, which is directly based on the kinematics of the fragmentation code used by SOPHIA (cf. Section 4.4.2).

3.3. Cross Sections

We parameterize the cross sections of photohadronic interactions following Mücke et al. (2000a). We split the processes into three parts: resonant, direct, and multi-pion production. The different contributions are shown in Figure 1.

The low-energy part of the total cross section is dominated by excitations and decays of baryon resonances N and Δ. The cross sections for these processes can be described by the Breit–Wigner formula:

$$\begin{eqnarray} \sigma _{\mathrm{BW}}^{\mathrm{IT}}(\epsilon _r) &=& \frac{s}{\big(s-m_p^2\big)^2}\frac{4\pi (2J+1)B_\gamma B_{\mathrm{out}}s\Gamma ^2}{(s-M^2)^2+s\Gamma ^2}\nonumber\\ &=& B_{\mathrm{out}}^{\mathrm{IT}} \frac{s }{\epsilon _r^2}\frac{\sigma _0^{\mathrm{IT}}(\Gamma ^{\mathrm{IT}})^2 s}{(s-(M^{\mathrm{IT}})^2)^2+(\Gamma ^{\mathrm{IT}})^2 s}, \end{eqnarray} \tag{ 20 }$$

with s(_r) given by Equation (5). Here, J, M, and Γ are the spin, the nominal mass, and the width of the resonance, respectively. We consider the energy-dependent branching ratios B_out and the resonances shown in Appendix B of Mücke et al. (2000a); the total branching ratios are also listed in Table 2. Note that the energy-dependent branching ratios respect the different energy thresholds for different channels (e.g., interaction types 1–4). For simplicity, we neglect the slight cross section differences between nγ and pγ interactions and use the values for protons, which implies that we do not take into account the N(1675) resonance.⁵ In Table 2, we list all of the resonances that were considered and their contributions to the different interaction types. To account for the phase-space reduction near threshold, the function Z^IT is introduced and multiplied in Equation (20):

$$\begin{equation} Z^{\mathrm{IT}}\big(\epsilon _r,\epsilon _{\mathrm{th}}^{\mathrm{IT}},w^{\mathrm{IT}}\big)= \left\lbrace \begin{array}{@{}ll@{}} 0& {{\rm if}\ \epsilon _r\le \epsilon _{\mathrm{th}}^{\mathrm{IT}},} \\ \frac{\epsilon _r-\epsilon _{\mathrm{th}}^{\mathrm{IT}}}{w^{\mathrm{IT}}}& {{\rm if}\ \epsilon _{\mathrm{th}}^{\mathrm{IT}}<\epsilon _r < w^{\mathrm{IT}}+\epsilon _{\mathrm{th}}^{\mathrm{IT}},}\\ 1& {{\rm if}\ \epsilon _r \ge w^{\mathrm{IT}}+\epsilon _{\mathrm{th}}^{\mathrm{IT}},} \end{array}\right. \end{equation} \tag{ 21 }$$

with ^IT_th and w^IT taking the values shown in Table 2.

Table 2. Summary of the Resonances Considered

Resonance	Properties			Z		Total Bout
	M (GeV)	Γ (GeV)	σ0 (μbarn)	th (GeV)	w (GeV)	IT 1	IT 2	IT 3
ΔIT(1232)	1.231	0.11	31.125	0.152	0.17	100%	⋅⋅⋅	⋅⋅⋅
NIT(1440)	1.440	0.35	1.389	0.152	0.38	67%	33%	⋅⋅⋅
NIT(1520)	1.515	0.11	25.567	0.152	0.38	52%	42%	6%
NIT(1535)	1.525	0.10	6.948	0.152	0.38	45%	⋅⋅⋅	⋅⋅⋅
NIT(1650)	1.675	0.16	2.779	0.152	0.38	75%	14%	⋅⋅⋅
NIT(1680)	1.680	0.125	17.508	0.152	0.38	64%	22%	14%
ΔIT(1700)	1.690	0.29	11.116	0.152	0.38	14%	55%	31%
ΔIT(1905)	1.895	0.35	1.667	0.152	0.38	14%	13%	73%
ΔIT(1950)	1.950	0.30	11.116	0.152	0.38	37%	39%	24%

Notes. Properties of the resonances and parameters for the quenching function (Equation (21)) taken from Mücke et al. (2000a). The four rightmost columns refer to possible interaction types IT in Table 1 and show the total branching ratios. In fact, at the end, each combination of a specific resonance and a particular decay chain corresponds to a separate interaction type, such as Δ₁(1232).

Download table as:  ASCIITypeset image

For the direct and the multi-pion cross section, we adopt the formulae from SOPHIA. The direct part is given by

$$\begin{eqnarray} \sigma ^{\mathrm{T}_1}(\epsilon _r)&=& 92.7 \, \mathrm{Pl}(\epsilon _r,0.152,0.25,2)+40\exp \left(\! -\frac{(\epsilon _r-0.29)^2}{0.002} \!\right)\nonumber\\ &&-15\exp \left(-\frac{(\epsilon _r-0.37)^2}{0.002} \right) \end{eqnarray} \tag{ 22 }$$

$$\begin{equation} \sigma ^{\mathrm{T}_2}(\epsilon _r)=37.7 \, \mathrm{Pl}(\epsilon _r,0.4,0.6,2) \end{equation} \tag{ 23 }$$

with

$$\begin{equation} \mathrm{Pl}(\epsilon _r,\epsilon _{\mathrm{th}},\epsilon _{\rm max},\alpha)= \left\lbrace \begin{array}{@{}ll@{}} 0& {{\rm if}\ \epsilon _r\le \epsilon _{\mathrm{th}}}\\ \ms \big(\frac{\epsilon _r-\epsilon _{\mathrm{th}}}{\epsilon _{\rm max}-\epsilon _{\mathrm{th}}}\big)^{A-\alpha }\big(\frac{\epsilon _r}{\epsilon _{\rm max}}\big)^{-A} &{\rm else,} \end{array}\right. \end{equation} \tag{ 24 }$$

where A = α_max/_th. The cross sections are given in μbarn, and the interaction types are taken from Table 1, i.e., T₁ is the single-pion direct production and T₂ is the two-pion direct production. The multi-pion cross section can be parameterized by the sum of the following two formulae (cross sections in μbarn):

$$\begin{equation} \sigma ^{\mathrm{Multi\hbox{-}\pi _1}}(\epsilon _r)= 80.3 \, Q_f(\epsilon _r,0.5,0.1) \, s^{-0.34}, \end{equation} \tag{ 25 }$$

$$\begin{eqnarray} \sigma ^{\mathrm{Multi\hbox{-}\pi _2}}(\epsilon _r)&=&\left[1-\exp \left(-\frac{\epsilon _r-0.85}{0.69}\right)\right]\nonumber\\ &&\times (29.3\, s^{-0.34}+59.3 \, s^{0.095})\ \ {\rm if\ \epsilon _r>0.85.}\qquad \end{eqnarray} \tag{ 26 }$$

We have compared our resonant and direct production multiplicity time cross sections and proton-to-neutron ratios with SOPHIA (cf. Figures 9–11 in Mücke et al. 2000a), and they are in excellent agreement.

Here we discuss simplified models based on Section 3, where we demand the features given in the introduction. For the resonances, we propose two methods; one is supposed to be more accurate and the other one is supposed to be faster.

4.1. Factorized Response Function

First of all, it turns out to be useful to simplify the response function in Equation (9). In our simplified approaches, we follow the following general idea: in Equations (7) and (9), in principle, (at least) two integrals have to be evaluated. Let us now split up the interactions in interaction types such that χ^IT(_r) ≡ χ^IT and M^IT_b(_r) ≡ M^IT_b are approximately constants independent of _r, but dependent on the interaction type. The response function in Equation (9) then factorizes in

$$\begin{eqnarray} R^{\mathrm{IT}}(x,y)&=& \delta (x - \chi ^{\mathrm{IT}}) \, M_b^{\mathrm{IT}} \, f^{\mathrm{IT}}(y) \quad {\rm with} \nonumber\\ f^{\mathrm{IT}}(y) &\equiv& \frac{1}{2y^2} \int _{\epsilon _{\mathrm{th}}}^{2y} d \epsilon _r \, \epsilon _r \, \sigma ^{\mathrm{IT}}(\epsilon _r) \,. \end{eqnarray} \tag{ 27 }$$

The part δ(x − χ^IT) describes at what energy the secondary particle is found, whereas the part M^IT_b f^IT(y) describes the production rate of the specific species b as a function of y. If only the number of produced secondary particles is important, it is often useful to show the effective F_b(y) ≡ ∑_ITM^IT_b f^IT(y).

As the next steps, we evaluate Equation (7) with Equation (27) by integrating over dx/x = −dE_p/E_p and by rewriting the ε-integral as y = E_pε/m_p integral:

$$\begin{equation} Q_b^{\mathrm{IT}} = N_p \left(\frac{E_b}{\chi ^{\mathrm{IT}}} \right) \frac{m_p}{E_b} \int _{\epsilon _{\mathrm{th}}/2}^\infty dy \, n_\gamma \left(\frac{m_p \, y \, \chi ^{\mathrm{IT}}}{E_b} \right) M_b^{\mathrm{IT}} \, f^{\mathrm{IT}}(y) \,. \end{equation} \tag{ 28 }$$

This (single) integral is relatively simple and fast to compute if the simplified response function f^IT(y) together with χ^IT is given. Therefore, the original double integral simplifies in this single integral, summed over a number of appropriate interaction types. In the following, we therefore define M^IT_b, f^IT(y), and χ^IT for suitable interaction types.

4.2. Resonances

Here we describe two alternatives for including the resonances. Model A is more accurate but slower to evaluate because it includes more interaction types, whereas model B is faster for computations since there are only two interaction types defined for the resonances.

4.2.1. Simplified Model A (Sim-A)

In this simplified resonance treatment, we make use of the fact that the resonances have cross sections peaked in _r. In principle, these can be approximated from Equations (20) and (21) by a δ-function

$$\begin{eqnarray} \sigma ^{\mathrm{IT}}(\epsilon _r) &\simeq & B^{\mathrm{IT}}_{\mathrm{out}} \, \hat{\Gamma }^{\mathrm{IT}} \, \delta \big(\epsilon _r-\epsilon _{r}^{\mathrm{IT}, 0}\big) \quad {\rm with}\nonumber\\ \quad \hat{\Gamma }^{\mathrm{IT}} &=& \int _0^{\infty } \sigma ^{\mathrm{IT}}_{\mathrm{BW}} \, Q_f^{\mathrm{IT}}\big(\epsilon _r,\epsilon _{\mathrm{th}}^{\mathrm{IT}},w^{\mathrm{IT}}\big) d \epsilon _r \,, \end{eqnarray} \tag{ 29 }$$

where $\hat{\Gamma }^{\mathrm{IT}}$ and ^IT,0_r correspond to surface area (in μbarn–GeV) and position (in GeV) of the resonance, and are process-dependent constants. Therefore, using Equation (29), the _r-integral in Equation (9) can be easily performed, and we obtain again the simplified Equation (27) with

$$\begin{eqnarray} f^{\mathrm{IT}}(y) &=& \frac{1}{2 y^2} \, \epsilon _{r}^{\mathrm{IT}, 0} \, B^{\mathrm{IT}}_{\mathrm{out}} \, \hat{\Gamma }^{\mathrm{IT}} \, \Theta \big(2 y - \epsilon _r^{\mathrm{IT},0}\big) \quad {\rm and}\nonumber\\ \quad \chi ^{\mathrm{IT}} &\equiv & \chi ^{\mathrm{IT}}\big(\epsilon _r^{\mathrm{IT}, 0}\big)\,. \end{eqnarray} \tag{ 30 }$$

The relevant parameters for all interaction types can be found in Table 3. Note that M^IT_π is the total (i.e., energy-independent) branching ratio from Table 2. In addition, note that in some cases, the resonance peak may be below the threshold for some of the interaction types. However, for reasonably broad energy distributions with which to be folded with, this should be a good approximation if the total branching ratios are used. The pion spectra are finally obtained from Equation (28), summing over all interactions listed in Table 3.

Table 3. Parameters for the Resonant Pion Production in Our Simplified Treatment A

IT	BITout	IT,0r	$\hat{\Gamma }^{\mathrm{IT}}$	χIT	Initial Proton			Initial Neutron
					$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$
Δ1(1232)	1.00	0.34	66.69	0.22	$\frac{2}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$
N1(1440)	0.67	0.64	4.26	0.29	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$
N2a(1440)	0.33	0.64	4.26	0.14	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{2}$	$\frac{1}{3}$	$\frac{1}{2}$	$\frac{1}{6}$
N2b(1440)	0.33	0.64	4.26	0.19	$\frac{1}{3}$	$\frac{11}{18}$	$\frac{1}{18}$	$\frac{1}{3}$	$\frac{1}{18}$	$\frac{11}{18}$
N1(1520)	0.52	0.75	27.01	0.31	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$
N2a(1520)	0.42	0.75	27.01	0.17	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{2}$	$\frac{1}{3}$	$\frac{1}{2}$	$\frac{1}{6}$
N2b(1520)	0.42	0.75	27.01	0.18	$\frac{1}{3}$	$\frac{11}{18}$	$\frac{1}{18}$	$\frac{1}{3}$	$\frac{1}{18}$	$\frac{11}{18}$
N3(1520)	0.06	0.75	27.01	0.22	$\frac{2}{3}$	1	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	1
N1(1535)	0.45	0.77	6.59	0.32	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$
N1(1650)	0.75	1.03	3.19	0.35	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$
N2a(1650)	0.14	1.03	3.19	0.23	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{2}$	$\frac{1}{3}$	$\frac{1}{2}$	$\frac{1}{6}$
N2b(1650)	0.14	1.03	3.19	0.17	$\frac{1}{3}$	$\frac{11}{18}$	$\frac{1}{18}$	$\frac{1}{3}$	$\frac{1}{18}$	$\frac{11}{18}$
N1(1680)	0.64	1.04	15.72	0.35	$\frac{1}{3}$	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$
N2a(1680)	0.22	1.04	15.72	0.24	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{2}$	$\frac{1}{3}$	$\frac{1}{2}$	$\frac{1}{6}$
N2b(1680)	0.22	1.04	15.72	0.17	$\frac{1}{3}$	$\frac{11}{18}$	$\frac{1}{18}$	$\frac{1}{3}$	$\frac{1}{18}$	$\frac{11}{18}$
N3(1680)	0.14	1.04	15.72	0.23	$\frac{2}{3}$	1	$\frac{1}{3}$	$\frac{2}{3}$	$\frac{1}{3}$	1
Δ1(1700)	0.14	1.05	21.68	0.35	$\frac{2}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$
Δ2a(1700)	0.55	1.05	21.68	0.24	$\frac{1}{15}$	$\frac{8}{15}$	$\frac{2}{5}$	$\frac{1}{15}$	$\frac{2}{5}$	$\frac{8}{15}$
Δ2b(1700)	0.55	1.05	21.68	0.16	$\frac{2}{5}$	$\frac{19}{45}$	$\frac{8}{45}$	$\frac{2}{5}$	$\frac{8}{45}$	$\frac{19}{45}$
Δ3(1700)	0.31	1.05	21.68	0.23	$\frac{1}{3}$	1	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	1
Δ1(1905)	0.14	1.45	3.03	0.38	$\frac{2}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$
Δ2a(1905)	0.13	1.45	3.03	0.29	$\frac{1}{15}$	$\frac{8}{15}$	$\frac{2}{5}$	$\frac{1}{15}$	$\frac{2}{5}$	$\frac{8}{15}$
Δ2b(1905)	0.13	1.45	3.03	0.15	$\frac{2}{5}$	$\frac{19}{45}$	$\frac{8}{45}$	$\frac{2}{5}$	$\frac{8}{45}$	$\frac{19}{45}$
Δ3(1905)	0.73	1.45	3.03	0.23	$\frac{1}{3}$	1	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	1
Δ1(1950)	0.37	1.56	16.46	0.39	$\frac{2}{3}$	$\frac{1}{3}$	⋅⋅⋅	$\frac{2}{3}$	⋅⋅⋅	$\frac{1}{3}$
Δ2a(1950)	0.39	1.56	16.46	0.30	$\frac{1}{15}$	$\frac{8}{15}$	$\frac{2}{5}$	$\frac{1}{15}$	$\frac{2}{5}$	$\frac{8}{15}$
Δ2b(1950)	0.39	1.56	16.46	0.15	$\frac{2}{5}$	$\frac{19}{45}$	$\frac{8}{45}$	$\frac{2}{5}$	$\frac{8}{45}$	$\frac{19}{45}$
Δ3(1950)	0.24	1.56	16.46	0.23	$\frac{1}{3}$	1	$\frac{2}{3}$	$\frac{1}{3}$	$\frac{2}{3}$	1

Notes. The units of ^IT,0_r are GeV and the units of $\hat{\Gamma }^{\mathrm{IT}}$ are μbarn GeV.

Download table as:  ASCIITypeset image

4.2.2. Simplified Model B (Sim-B)

In this model, we take into account the width of the resonances and approximate them by constant cross sections within certain energy ranges. Compared to approaches such as in Reynoso & Romero (2009), we distinguish between π⁺ and π⁻ production (i.e., do not add these fluxes) and include the effects from the higher resonances.

From Figure 1, one can read off that there are two classes of resonances: The first peak in Figure 1 is dominated by the Δ(1232)-resonance (lower resonance—LR), whereas the higher resonances (HR) contribute at larger energies. In addition, the kinematics (cf. χ in Table 3) and the multiplicities (cf. Table 1) are very different. For example, protons and neutrons interactions via the Δ(1232)-resonance happen through only one interaction type, and produce either only π⁺ or only π⁻ (and π⁰). For the higher resonances, the pions are produced in all pion charges.

For the resonances, we define two interaction types LR and HR, as shown in Figure 2 (boxes). The interaction type LR corresponds to Δ₁(1232), whereas the interaction type HR contains the higher resonances. The properties of these interaction types are summarized in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Parameters for the Resonant Pion Production in Our Simplified Treatment B

IT	r Range (GeV)	σ (μbarn)	Initial Proton				Initial Neutron
			$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	Mn	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	Mn	K
			$\chi _{\pi ^0}$	$\chi _{\pi ^+}$	$\chi _{\pi ^-}$	Mp	$\chi _{\pi ^0}$	$\chi _{\pi ^+}$	$\chi _{\pi ^-}$	Mp
LR	0.2–0.5	200	2/3	1/3	⋅⋅⋅	1/3	2/3	⋅⋅⋅	1/3	2/3	0.22
			0.22	0.22	⋅⋅⋅	2/3	0.22	⋅⋅⋅	0.22	1/3
HR	0.5–1.2	90	0.47	0.77	0.34	0.43	0.47	0.34	0.77	0.57	0.39
			0.26	0.25	0.22	0.57	0.26	0.22	0.25	0.43

Download table as:  ASCIITypeset image

The surface area covered by these cross sections, corresponding to the $\hat{\Gamma }$ in simplified model A, is chosen for LR and HR such that the pion spectra match to these of the previous section for typical power-law spectra (with spectral indices of about two). Of course, there can never be an exact match because the contributions from the individual resonances depend on the spectral index. However, as we will demonstrate, this estimate is good enough for our purposes.

The averaged numbers for the multiplicities and inelasticities for the higher resonance interaction type are estimated from the interaction rate Equation (4) by assuming that all resonances contribute simultaneously and by weighting the interaction type-dependent part $\epsilon _{r}^{\mathrm{IT}, 0} \, \hat{\Gamma }^{\mathrm{IT}} \, B^{\mathrm{IT}}_{\mathrm{out}}$ (cf. Equation (B3) using Equation (30); cf. Appendix B). By the same procedure we obtain the average χ-values from Equation (9). Note that the π⁻ (for proton interactions) are, in average, reconstructed at lower energies because these are mostly produced by interaction type 2, which has smaller χ-values (cf. Table 3). In addition, note that the total number of charged pions is, in fact, close to one per interaction here, and the total number of pions is close to 1.5, since in interaction types 2 and 3 more than one pion is produced.

The pion spectra are then computed from Equation (28) with

$$\begin{eqnarray} &&f^{\mathrm{LR}} =\nonumber\\ \ms &&\ \ \ \left\lbrace \begin{array}{@{}ll@{}} 0 & 2y < 0.2 \, \mathrm{GeV} \\ 200 \, \mu \mathrm{barn} \, \big(1- \frac{(0.2 \, \mathrm{GeV})^2}{(2y)^2} \big)& 0.2 \, \mathrm{GeV} \le 2y <0.5 \, \mathrm{GeV}\\ \ms 200 \, \mu \mathrm{barn} \, \big(\frac{(0.5 \, \mathrm{GeV})^2-(0.2 \, \mathrm{GeV})^2}{(2y)^2} \big) & 2y \ge 0.5 \, \mathrm{GeV}\\ \end{array} \right.\nonumber\\ \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&f^{\mathrm{HR}} =\nonumber\\ &&\ \ \ \left\lbrace \begin{array}{@{}ll@{}} 0 & 2y < 0.5 \, \mathrm{GeV} \\ 90 \, \mu \mathrm{barn} \, \big(1- \frac{(0.5 \, \mathrm{GeV})^2}{(2y)^2} \big) & 0.5 \, \mathrm{GeV} \le 2y < 1.2 \, \mathrm{GeV}\\ \ms 90 \, \mu \mathrm{barn} \, \big(\frac{(1.2 \, \mathrm{GeV})^2-(0.5 \, \mathrm{GeV})^2}{(2y)^2} \big) & 2y \ge 1.2 \, \mathrm{GeV} \,, \end{array} \right.\nonumber\\ \end{eqnarray} \tag{ 32 }$$

as can be shown from Equation (27). Here, M^IT_π and χ^IT (needed in Equation (28)) are chosen according to interaction type, initial proton or neutron, and final pion, as given in Table 4.

4.2.3. Comparison of Resonance Response Functions

In Figure 3, we compare the response function F_π(y) ≡ ∑_ITM^IT_π f^IT(y) between simplified model Sim-A (thin solid curves) and Sim-B (thick curves), summed over the resonances. This function is proportional to the number of produced pions of a certain species as a function of y, whereas the x-dependent part in Equation (27) describes the energy at which the pions are found. Obviously, the function is much smoother for Sim-B than for Sim-A. Because of only a few contributing interaction types and the smoothness of the function, the evaluation will be much faster. Once the photon spectrum is folded in, the contributions of both response functions will be very similar.

Zoom In Zoom Out Reset image size

Note that model Sim-B includes the part of the Δ(1232)-resonance below the peak. This can be seen by comparison with the dashed curve, which represents model Sim-A but the interaction type Δ₁(1232) replaced by the full Breit–Wigner form.

4.3. Direct Production

For direct production, we also follow the approach in Section 4.1. However, this interaction type is tricky, since the kinematic function χ^IT(_r) is strongly dependent on the interaction energy (see Appendix A for details). This implies that the δ-distribution in Equation (6) is not a good approximation and needs to be replaced by a broader distribution function. In addition, the distribution of scattering angles will lead to smearing effects. In this case, it turns out to be a good approximation for splitting the direct production in different interaction types with different characteristic values of χ^IT as a function of _r, which simulates such a broad distribution after the integration of the input spectra. We define three interaction types T_1L, T_1M, and T_1H for direct one-pion production and four for two-pion production, namely T_2aL, T_2aM, and T_2aH for the first pion, and T_2b for the second pion. The interaction types are shown in Table 5, and the names correspond to Table 1, split into low (L), medium (M), and high (H) energy parts in _r, respectively, limited by ^IT_min and ^IT_max. After this splitting, the additional effect of the scattering angle smearing turns out to be small, as we will demonstrate later. We obtain the function f^IT for these interaction types from Equation (27) as

$$\begin{equation} f^{\mathrm{IT}} = \left\lbrace \begin{array}{@{}ll@{}} 0 & 2y < \epsilon _{\mathrm{min}}^{\mathrm{IT}} \\ \ms \frac{1}{2y^2}\big(I^{\mathrm{IT}}(2y)-I^{\mathrm{IT}}\big(\epsilon _{\mathrm{min}}^{\mathrm{IT}}\big)\big) & \epsilon _{\mathrm{min}}^{\mathrm{IT}} \le 2y < \epsilon _{\mathrm{max}}^{\mathrm{IT}}\\ \ms \frac{1}{2y^2}\big(I^{\mathrm{IT}}\big(\epsilon _{\mathrm{max}}^{\mathrm{IT}}\big)-I^{\mathrm{IT}}\big(\epsilon _{\mathrm{min}}^{\mathrm{IT}}\big)\big) & 2y \ge \epsilon _{\mathrm{max}}^{\mathrm{IT}},\\ \end{array} \right.\\ \end{equation} \tag{ 33 }$$

with I^IT(2y) re-parameterized using x ≡ log₁₀(y/GeV):

$$\begin{eqnarray} I^{\mathrm{T}_{1}}(2y) & = & \left\lbrace \begin{array}{@{}l@{}} 0 \qquad\qquad\quad \ \ \ 2 y < 0.17 \, \mathrm{GeV} \\ \ms 35.9533+84.0859 x+110.765 x^2+102.728 x^3\\ \quad+\,40.4699 x^4 \quad \ \ 0.17\, \mathrm{GeV}< 2 y<0.96 \, \mathrm{GeV} \\ \ms 30.2004 + 40.5478 x + 2.03074 x^2 - 0.387884 x^3\\ \quad +\, 0.025044 x^4 \quad 2y > 0.96 \, \mathrm{GeV} \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} I^{\mathrm{T}_{2}}(2y) & = & \left\lbrace \begin{array}{@{}l@{}} 0 \qquad\qquad\qquad\qquad 2 y < 0.4 \, \mathrm{GeV} \\ \ms -3.4083+\frac{16.2864}{2\,y}\\ \quad +\,40.7160 \ln {(2 y)} \qquad \ 2 y \ge 0.4 \, \mathrm{GeV.} \\ \end{array} \right.\qquad \end{eqnarray} \tag{ 35 }$$

Note that ^IT_min, ^IT_max, the multiplicities, and χ^IT for the different interaction types can be read off from Table 5. The integral value I^IT is the same for interaction types T_1L, T_1M, T_1H (Equation (34)), and for interaction types T_2aL, T_2aM, T_2aH, and T_2b (Equation (35)). In Equation (28), all the interaction types in Table 5 have to be summed.

Table 5. Parameters for the Direct Pion Production in Our Simplified Treatment

IT	ITmin (GeV)	ITmax (GeV)	χ	K	Initial Proton			Initial Neutron
					$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$
T1L	0.17	0.56	0.13	0.13	⋅⋅⋅	1	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	1
T1M	0.56	10	0.05	0.05	⋅⋅⋅	1	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	1
T1H	10	∞	0.001	0.001	⋅⋅⋅	1	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	1
T2aL	0.4	1.58	0.08	0.28	⋅⋅⋅	$\frac{1}{4}$	$\frac{3}{4}$	⋅⋅⋅	$\frac{3}{4}$	$\frac{1}{4}$
T2aM	1.58	10	0.02	0.22	⋅⋅⋅	$\frac{1}{4}$	$\frac{3}{4}$	⋅⋅⋅	$\frac{3}{4}$	$\frac{1}{4}$
T2aH	10	∞	0.001	0.201	⋅⋅⋅	$\frac{1}{4}$	$\frac{3}{4}$	⋅⋅⋅	$\frac{3}{4}$	$\frac{1}{4}$
T2b	0.4	∞	0.2	⋅⋅⋅	$\frac{1}{6}$	$\frac{3}{4}$	$\frac{1}{12}$	$\frac{1}{6}$	$\frac{1}{12}$	$\frac{3}{4}$

Notes. The inelasticity for T_2b is included in T_2a.

Download table as:  ASCIITypeset image

4.4. Multi-pion Production

Here we show two different approaches to the multi-pion channel. The first simplified approach will later be shown as model Sim-C; the second more refined approach will be used in all other models.

4.4.1. Simplified Kinematics

We start with the simplest example for multi-pion production, for which we follow Section 3.2. We assume χ^Multi-π ≃ 0.2, i.e., the response function factorizes as R^Multi-π = δ(x − 0.2) M^Multi-π_π f^Multi-π(y). We obtain from Equations (27) and (28)

$$\begin{eqnarray} Q_{\pi }^{\mathrm{Multi\hbox{-}\pi }}(E_\pi) &=& N_p(5 E_\pi) \, \frac{m_p}{E_\pi } \, \int _{\epsilon _{\mathrm{th}}/2}^\infty d y \, M_\pi ^\mathrm{Multi\hbox{-}\pi } \nonumber\\ &&\times\! f^{\mathrm{Multi\hbox{-}\pi }}(y) \, n_\gamma \left(\frac{m_p y}{5 E_\pi } \right), \end{eqnarray} \tag{ 36 }$$

with

$$\begin{equation} f^{\mathrm{Multi\hbox{-}\pi }}(y) \equiv \frac{1}{2 y^2} \int _{\epsilon _{\mathrm{th}}}^{2y} d \epsilon _r \, \epsilon _r\, \sigma ^{\mathrm{Multi\hbox{-}\pi }}(\epsilon _r) \,. \end{equation} \tag{ 37 }$$

As described in Section 3, we assume that $M_{\pi ^0}=1$ ($M_{\pi ^0}=1$), $M_{\pi ^+}=1.2$ ($M_{\pi ^+}=0.85$), and $M_{\pi ^-} = 0.85$ ($M_{\pi ^-} =1.2$) for initial protons (neutrons). The function f^Multi-π(y) can be obtained using the sum of Equations (25) and (26). We numerically integrate it and re-parameterize it with x = log₁₀(y/GeV) by

$$\begin{equation} \frac{\hat{f}(x)_{\pi ^{+,0}}}{\mu {\rm barn}} = \left\lbrace \begin{array}{@{}l@{}} 0 \qquad\qquad\qquad\quad \, 2 y<0.5 \, \mathrm{GeV} \\ \ms 87.5538 + 120.894 x - 98.4187 x^2 - 59.6965 x^3 \\ \quad+\, 67.2251 x^4 \qquad\ \ \ \, 0.5 \, \mathrm{GeV} \le 2 y \le 14 \, \mathrm{GeV} \\ \ms 131.839 - 25.3296 x + 10.612 x^2 - 0.858307 x^3 \\ \quad+\, 0.0493614 x^4 \qquad 2y > 14 \, \mathrm{GeV} \end{array} \right.\nonumber\\ \end{equation} \tag{ 38 }$$

for π⁰, and π⁺ (initial proton) or π⁻ (initial neutron), and

$$\begin{eqnarray} &&\frac{\hat{f}(x)_{\pi ^{-}}}{\mu {\rm barn}} =\nonumber\\ &&\quad \left\lbrace \begin{array}{@{}l@{}} 0 \qquad\qquad\qquad\qquad 2 y< 1 \, \mathrm{GeV} \\ \ms 73.9037 + 187.526 x - 161.587 x^2 - 206.268 x^3 \\ \quad +\, 354.02 x^4 - 129.759 x^5 \quad 1 \, \mathrm{GeV} \le 2 y \le 10 \, \mathrm{GeV} \\ \ms 131.839 - 25.3296 x + 10.612 x^2 - 0.858307 x^3 \nonumber\\ \quad+\, 0.0493614 x^4 \qquad 2y > 10 \, \mathrm{GeV} \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 39 }$$

for π⁻ (initial proton) or π⁺ (initial neutron). Note that the different function for π⁻ comes from the different threshold below which we have set the cross section to zero, as discussed in Section 3.2. For y ≫ 10⁴ GeV, this function still increases as extrapolation of Figure 1 (the cross section is not measured in that energy range).

4.4.2. Kinematics Simulated from SOPHIA

Similar to direct production, the δ-distribution in Equation (6) is not a good approximation and needs to be replaced by a broader distribution function. We again solve this by splitting the multi-pion production in different interaction types with different characteristic values of χ^IT, which simulates such a broad distribution after the integration of the input spectra. Compared to direct production, the cross section can be parameterized relatively easily by step functions. However, the pion multiplicities in fact increase with energy. This means that the higher the interaction energy is, the more pions are produced, which (in average) are found at lower energies. In addition, the final energy distribution functions are broad.

We solve this strongly scale-dependent behavior by dividing the _r-range into seven interaction types M_i, each with a particular average cross section and average pion multiplicities, which we take directly from SOPHIA. In addition, we split the pions in each of these samples into lower energy (L) and higher energy (H) parts, which are reconstructed at different values of χ^IT to simulate the breadth of the distributions within the same _r. Typically, the L-sample corresponds to the peak of the distribution (at least for high energies), whereas the H-sample simulates the tail of the distribution. The splitting of the multiplicities can be performed automatically once a splitting point is defined. We choose the χ-values of the interaction types to simulate the peaks of the distribution and to reproduce the total energy going into pions in SOPHIA. This treatment is, of course, not extremely accurate, but it allows us to use the fast approach in Equation (27) while obtaining accuracy at high energies.

The function f^M can be easily calculated from Equation (27) since the cross section is assumed to be constant within each IT:

$$\begin{equation} f^{\mathrm{M}_i} = \left\lbrace \begin{array}{@{}ll@{}} 0 & 2y < \epsilon _{\mathrm{min}}^{\mathrm{M}_i} \\ \ms \frac{\sigma ^{\mathrm{M}_i}}{(2y)^2}\big((2y)^2-\big(\epsilon _{\mathrm{min}}^{\mathrm{M}_i}\big)^2 \big) & \epsilon _{\mathrm{min}}^{\mathrm{M}_i} \le 2y < \epsilon _{\mathrm{max}}^{\mathrm{M}_i}\\ \ms \frac{\sigma ^{\mathrm{M}_i}}{(2y)^2}\big(\big(\epsilon _{\mathrm{max}}^{\mathrm{M}_i}\big)^2-\big(\epsilon _{\mathrm{min}}^{\mathrm{M}_i}\big)^2 \big) & 2y \ge \epsilon _{\mathrm{max}}^{\mathrm{M}_i}.\\ \end{array} \right.\\ \end{equation} \tag{ 40 }$$

The interaction types are listed in Table 6, where the parameters for this equation can be found, as well as the χ and K values and multiplicities for the individual interaction types.

Table 6. Parameters for Multi-pion Production in Our SOPHIA-based Treatment

IT	ITmin (GeV)	ITmax (GeV)	σ	χ	K	Initial Proton			Initial Neutron
						$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$	$M_{\pi ^0}$	$M_{\pi ^+}$	$M_{\pi ^-}$
M1L	0.5	0.9	60	0.1	0.27	0.32	0.34	0.04	0.32	0.04	0.34
M1H	0.5	0.9	60	0.4	⋅⋅⋅	0.17	0.29	0.05	0.17	0.05	0.29
M2L	0.9	1.5	85	0.15	0.34	0.42	0.31	0.07	0.42	0.07	0.31
M2H	0.9	1.5	85	0.35	⋅⋅⋅	0.19	0.35	0.08	0.19	0.08	0.35
M3L	1.5	5.0	120	0.15	0.39	0.59	0.57	0.30	0.59	0.30	0.57
M3H	1.5	5.0	120	0.35	⋅⋅⋅	0.16	0.21	0.13	0.16	0.13	0.21
M4L	5.0	50	120	0.07	0.49	1.38	1.37	1.11	1.38	1.11	1.37
M4H	5.0	50	120	0.35	⋅⋅⋅	0.16	0.25	0.23	0.16	0.23	0.25
M5L	50	500	120	0.02	0.45	3.01	2.86	2.64	3.01	2.64	2.86
M5H	50	500	120	0.5	⋅⋅⋅	0.20	0.21	0.14	0.20	0.14	0.21
M6L	500	5000	120	0.007	0.44	5.13	4.68	4.57	5.13	4.57	4.68
M6H	500	5000	120	0.5	⋅⋅⋅	0.27	0.29	0.12	0.27	0.12	0.29
M7L	5000	∞	120	0.002	0.44	7.59	6.80	6.65	7.59	6.65	6.80
M7H	5000	∞	120	0.6	⋅⋅⋅	0.26	0.27	0.13	0.26	0.13	0.27

Notes. The inelasticities for the H-sample are included in the L-sample.

Download table as:  ASCIITypeset image

4.5. Kaon Production

Here we include K⁺ production in our simplified model. If needed, other (sub-leading) channels can be added to our framework as well (such as K⁻ production if the neutrino–antineutrino ratio in the higher energy regime is studied). This example should also serve as an illustration for how to add processes to our model. Note that this implementation, however, is more primitive than the pion production above.

The production of K⁺ in photohadronic interactions cannot be assigned to a single resonance, but comes from a number of resonances mainly decaying into K⁺ and Λ or Σ⁰ (with relatively similar masses, which means that we can neglect the mass difference). In addition, at high interaction energies, multi-fragmentation, similar to multi-pion production, contributes significantly. The production cross section up to _r ≃ 2 GeV has been measured (Saphir Collaboration 1998). For higher energies, where no data are available, one may use extrapolations from models (Lee et al. 2001; Asano & Nagataki 2006). Approximating the data in Lee et al. (2001) and extrapolating according to Lee et al. (2001), we define an interaction type "KP" and model the total K⁺ production cross section as (cf. dashed curve in Figure 2):

$$\begin{eqnarray} \sigma ^{\mathrm{KP}} & = & \left\lbrace \begin{array}{@{}ll@{}} 0 & \epsilon _r < 1.0 \, \mathrm{GeV} \\ \ms 2.0 \, \mu \mathrm{barn} & 1.0 \, \mathrm{GeV}\le \epsilon _r < 1.2 \, \mathrm{GeV}\\ \ms 3.7 \, \mu \mathrm{barn} & 1.2 \, \mathrm{GeV}\le \epsilon _r < 1.65 \, \mathrm{GeV}\\ \ms 2.7 \, \mu \mathrm{barn} & \epsilon _r \ge 1.65 \, \mathrm{GeV}. \\ \end{array} \right. \end{eqnarray} \tag{ 41 }$$

Accordingly, we obtain the response function from Equation (27) with $M_{K^+}=1$:

$$\begin{eqnarray} &&f^{\mathrm{KP}} =\nonumber\\ &&\ \ \ \left\lbrace \begin{array}{ll} 0 & 2y < 1.0 \, \mathrm{GeV} \\ \ms 2.0 \, \mu \mathrm{barn} \, \big(1- \frac{1.0 \, \mathrm{GeV}^2}{(2y)^2} \big) & 1.0 \, \mathrm{GeV} \le 2y < 1.2 \, \mathrm{GeV}\\ \ms 3.7 \, \mu \mathrm{barn} \, \big(1- \frac{0.3 \, \mathrm{GeV}^2}{(2y)^2} \big) & 1.2 \, \mathrm{GeV} \le 2y < 1.65 \, \mathrm{GeV}\\ \ms 2.7 \, \mu \mathrm{barn} \, \big(1- \frac{0.16 \, \mathrm{GeV}^2}{(2y)^2} \big) & 2y \ge 1.65 \, \mathrm{GeV}.\\ \end{array} \right.\nonumber\\ \end{eqnarray} \tag{ 42 }$$

From Section 3.2, we have $\chi _{K^+} \simeq 0.35$ (computed for _r ≃ 1.4 GeV at the peak). The effects on primary cooling and secondary re-injection in Appendix B are negligible. Note that the absolute normalization of the kaons at high energies crucially depends on the extrapolation of the cross sections to high energies. Our kaon production is about a factor of two smaller at higher energies than in Lipari et al. (2007), since we assume the cross section to be constant at high energies.

4.6. Comparison of Individual Contributions

In Figure 4, we show the individual contributions to F_π(y) ≡ ∑_ITM^IT_π f^IT(y) for pion production as a function of y, which is proportional to the maximal available center-of-mass energy. This quantity describes the total number of pions of a particular species produced as a function of y, including multiplicity and cross section. However, it does not include the energy where the pions are found. For initial protons, π⁺ are most abundantly produced because of the direct production contribution, followed by π⁰ and then π⁻. For π⁺ production close to the threshold, the direct production and resonances are most important, whereas for π⁰ production, the direct production hardly has any impact, and the resonances dominate. For π⁻ production, at the threshold all processes contribute almost equally, only the LR does not take part. For larger y, multi-pion production dominates in all cases.

Zoom In Zoom Out Reset image size

As far as the relative contribution is concerned, thanks to direct production, the π⁺ is always produced most abundantly. As one can see from Figure 4, this statement is independent of the photon or proton spectra used because the response function for π⁺ is always larger than the one for π⁰ and in a large part of the energy range it is even about 50% larger. This means that for arbitrary proton and photon spectra, thanks to direct production, the π⁺:π⁰ ratio lies between about 1:1 and 3:2, whereas the Δ(1232)-approximation in Equation (1) predicts 1:2. In fact, one can show that the minimum of the charged to neutral pion ratio with respect to y is

$$\begin{equation} {\mathrm{min}_{y}} \frac{F_{\pi ^+}(y)+F_{\pi ^-}(y)}{F_{\pi ^0}(y)} \simeq 1.2 \, \end{equation} \tag{ 43 }$$

for arbitrary input spectra; see Mücke et al. (2000b) for a specific photon field.⁶ From this equation, any neutrino flux computed with the Δ(1232)-approximation and normalized to the observed photon spectrum, if mainly coming from π⁰ decays, is underestimated by a factor of at least 1.2/0.5 = 2.4, i.e., the neutrino flux should be about a factor of 2.4 larger. Even the often used approximation that 50% of all photohadronic interactions result in charged pions underestimates the neutrino flux by at least 20%. These numbers are to be interpreted as lower limits for the neutrino flux underestimation; the exact values depend on the input spectra and may be even higher.

Of course, the overall impact of the individual contributions depends on the proton and photon spectra with which the response function is to be folded. In order to check the impact of the individual contributions on typical AGN, GRB, or high-temperature blackbody (BB) spectra, we define three benchmark spectra in Appendix C. Note that all of the spectra shown in this work are given in the SRF. In Figure 5, we show the π⁺, π⁻, and π⁰ spectra for these benchmarks. The upper panels are for the GRB benchmark, the middle ones are for the AGN benchmark and the lower ones are for the BB benchmark. For the GRB benchmark, direct production dominates at low energies for π⁺ production, whereas π⁰ production at low energies is determined by the resonances. At all energies π⁻ production is dominated by multi-pion production, like the other spectra at high energies. Therefore, all different processes are important, but they contribute entirely differently to the different pion polarities. One can also see from this figure that the characteristic shape of the GRB pion or neutrino spectra often shown in the literature (resonance curves), which is flat in the middle energy range, is tilted upward due to multi-pion production. For the AGN benchmark (middle panels), π⁰ production is governed by the resonances in the whole energy range, π⁻ production is governed by the multi-pion production, and the π⁺ production is governed by similar contributions of all processes including direct production. In this case, the Δ-resonance approximation is more accurate in terms of the shapes of the spectra. However, other processes quantitatively contribute as well. For the BB benchmark (lower panels), π production is dominated by multi-pion production for most of the energies. Charged pion production is governed by direct and resonant processes only in the low-energy regions, and π⁰ production is governed by the resonances. In this case, ignoring the high-energy processes would lead to clearly misleading results.

Zoom In Zoom Out Reset image size

As far as the comparison among the different polarities is concerned, see Figure 6. For all spectra, the π⁺ are always most abundantly produced, as predicted above, followed by π⁰ and then π⁻. For the GRB benchmark at high energies, where multi-pion production dominates, the spectra are closest to each other. At low energies, there are significantly less π⁻ produced than the other two polarities. However, note that, thanks to multi-pion production, π⁻ is only suppressed by a factor of a few. The kaons, on the other hand, contribute about one to two orders less to the total meson fluxes. Nevertheless, there can be interesting effects in the high-energy regime if cooling effects are present. For the AGN benchmark because of the lower maximal proton energy times photon energy, even for large energies π⁻ is strongly suppressed. Otherwise, the result is qualitatively similar. For the BB spectrum we can nicely see the effect of lower multiplicities of π⁻ compared to π⁺ and π⁰ for high interaction energies (see Table 6) in the spectrum at energies higher than 10⁹ GeV. Otherwise, we have similar results as for the GRB benchmark.

Zoom In Zoom Out Reset image size

Here we compare the results of our simplified models with each other and with SOPHIA. First, we focus on the primaries produced in the photohadronic interactions, mostly the pions. Then we compute and compare the neutrino spectra. In all cases, we use initial protons for the photohadronic interactions.

5.1. Pion Spectra

The models considered in this subsection are listed in Table 7, where the (computational) complexity decreases from the top to the bottom. "SOPHIA" represents the output of the SOPHIA software, computed for our benchmark spectra from Appendix C. The computation with SOPHIA is described in Appendix D. The model "BW" corresponds to the basic physics of SOPHIA as described in Section 3, including double integration for direct production and multi-pion production from Section 4.4.2. In this case, we use Equations (7) and (9) for secondary production using the two-dimensional response function including all resonances with full Breit–Wigner forms and direct production as described in Appendix A. Compared to SOPHIA, the kinematics are considerably simplified. The multi-pion production is taken from Section 4.4.2 close to SOPHIA. "Sim-A" and "Sim-B" are the simplified models from Section 4, using the factorized response function introduced in Section 4.1. Whereas Sim-A treats all interaction types with the resonances explicitly, Sim-B defines an interaction type for the Δ(1232)-resonance, and one for the higher resonances. Therefore, Sim-B uses considerably less interaction types. Note that we have obtained Sim-A and Sim-B by condensing the information from BW stepwise. Both models use the multi-pion production from Section 4.4.2, whereas Sim-C is a simplified version of Sim-B including the multi-pion production in Section 4.4.1.

Table 7. Models Considered for the Comparison of Approaches

Abbrev.	Description	Complexity	Section/Refs.
SOPHIA	SOPHIA software with full kinematics	Monte Carlo method, ∼3000 × TSim-C	Mücke et al. (2000a)
BW	Resonances with full Breit–Wigner description; direct production with full response function; simple kinematics	Double integration, 45 IT, ∼120 × TSim-C	Section 3 incl. multi-pion from Section 4.4.2
Sim-A	Simplified model with factorized response function and resonance treatment A (δ-function approximation)	Single integration, 49 IT, ∼(3–4) × TSim-C	Section 4 with Section 4.2.1 and Section 4.4.2
Sim-B	Simplified model with factorized response function and resonance treatment B (step function approximation)	Single integration, 23 IT, ∼(2–3) × TSim-C	Section 4 with Section 4.2.2 and Section 4.4.2
Sim-C	Simplified model with factorized response function and resonance treatment B; simplified multi-pion production	Single integration, 10 IT, TSim-C	Section 4 with Section 4.2.2 and Section 4.4.1

Notes. The (computational) complexity decreases from the top to the bottom. The time needed for the computation of the photohadronics in model Sim-C is given by T_Sim-C. The comparison of the computation times is done for power-law spectra. For the computation with SOPHIA, we used 100,000 trials per proton bin.

Download table as:  ASCIITypeset image

We compare the pion spectra produced by these different models in Figure 7 (GRB benchmark), Figure 8 (AGN benchmark), and Figure 9 (blackbody (BB) benchmark). In these figures, the upper rows show the pion spectra for π⁺, π⁻, and π⁰ explicitly, whereas the lower rows show the pion ratios π⁺/π⁻, π⁺/π⁰, and π⁻/π⁰. Since the pion spectra for model BW, Sim-A, and Sim-B are so close to each other that they cannot be distinguished, we plot only the results of SOPHIA, model Sim-B, and model Sim-C in the upper rows of Figures 7–9. Whereas the charged pions lead to neutrino production, the neutral pions lead to photons. Therefore, the ratio π⁺/π⁰ determines, to leading order, the ratio between neutrinos and photons, which also often enters the computation of neutrino flux limits. On the other hand, the ratio π⁺/π⁻ affects the (electron) neutrino–antineutrino ratio. The ratio π⁻/π⁰ is shown for completeness. Note that the normalization of the different spectra is not chosen arbitrarily, but consistently to be able to directly compare the spectra; cf. Appendix D.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most importantly, as can be seen from the upper rows of the three figures, all the spectra of our simplified models (apart from Sim-C, may be) match the output from SOPHIA very well, both normalization and spectral shape. At high energies, however, where we imposed a sharp spectral cutoff, the spectra from SOPHIA are smeared out because of the more refined kinematics treatment, which can best be seen in Figure 9 for the BB benchmark where we have a sharp spectral cutoff in the proton spectrum. This difference is unavoidable, and the price that must be paid for an efficient simplified model. Nevertheless, note that in more realistic spectra, or spectra averaged over different sources, such sharp features in the spectra may not be present. Although it seems that there are hardly any differences between SOPHIA and our models at lower energies, one can read off from the pion ratio plots in the lower rows (on a linear vertical scale) that there are small deviations. Whereas the differences at very low and high energies, where the spectra rapidly break off, are not very surprising, there are some differences coming from the different kinematics treatment. For example, SOPHIA actually produces even about 20% more π⁺ than π⁰ in the lower energy range for the GRB benchmark and the middle energy range for the AGN benchmark (cf. middle lower panels). We have checked that this difference comes from neither the resonance nor direct production treatment; it also does not come from the relative contribution of both processes. Instead, for some processes, our kinematics are a bit oversimplified, since for these a considerable amount of pions is (in SOPHIA) reconstructed at lower energies. Another example for an obvious difference is the high-energy π⁺/π⁻ and π⁺/π⁰ difference in the lower panel of Figure 9 for the BB benchmark, the most challenging one for the treatment of multi-pion production. The discrepancy is a result of the simplified kinematic treatment of taking the same χ-values for π⁺, π⁰, and π⁻. This mainly effects the π⁻ spectra because they have, in comparison to π⁺ and π⁰, slightly different kinematics in SOPHIA. In the upper panel in Figure 9, a double-hump structure can be seen which follows from the kinematics treatment of multi-pion production that one part of the pions is reconstructed at lower energies (small χ) and the other at higher energies (large χ). If one averages over larger energy scales (such as in diffuse fluxes), such kinematic effects average out.

As far as the comparison among our simplified models is concerned, the differences are small compared to the effects of kinematics discussed above. In fact, model "BW," which was originally designed as the most accurate reproduction of SOPHIA, produces the results most different from SOPHIA, especially for the lower half of the energy range. The reason may be that the errors introduced by the approximations in Sim-A and Sim-B partly compensate for the errors from the simplified multi-pion production and direct channels. Model Sim-A was obtained from BW by assuming that all resonances are strongly peaked, whereas Sim-B was derived from this assumption by collecting the properties of the higher resonances into one interaction type. In the comparison of model Sim-B to Sim-C, the differences between the kinematics for multi-pion production from Section 4.4.2 and the simplified one from Section 4.4.1 can be seen nicely. For the GRB (see Figure 7) and the AGN benchmark (see Figure 8) this mainly affects the high-energy region, whereas in the BB case (see Figure 9) most of the energy range is affected. Since the computation time for Sim-B, for which the results do not differ significantly from model Sim-A, is only about a factor of two to three longer and the high-energy treatment is much more accurate than for Sim-C (especially close to the peaks), in the following section we focus on model Sim-B. As one can see in Figures 7–9, Sim-B is most accurate for power laws, which are our main interest. Compared to SOPHIA, we gain a factor of about 1000 (if implemented in C, depending on the integration method) in computation time for 100,000 trials per proton bin in SOPHIA (as we use for the GRB benchmark). The Sim-B spectra do not have any small wiggles because the computation is exact. However, note that the complexity of Sim-B increases with the number of interaction types, whereas the complexity of SOPHIA increases with the number of trials (and required smoothness of the functions).

5.2. Neutrino Spectra

In this section, we first review the production of neutrinos from pion and subsequent muon decays, as well as kaon decays. For the sake of completeness, we include the neutrinos from neutron decays, where the neutrons are produced by the photohadronic interactions. Then we compare our results to SOPHIA. We focus on model Sim B from the previous section only.

The neutrinos are mostly produced in the following two decay chains:

$$\begin{eqnarray} \pi ^+ & \rightarrow & \mu ^+ + \nu _\mu \,,\nonumber \\ & & \mu ^+ \rightarrow e^+ + \nu _e + \bar{\nu }_\mu \,, \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \pi ^- & \rightarrow & \mu ^- + \bar{\nu }_\mu \,, \nonumber \\ & & \mu ^- \rightarrow e^- + \bar{\nu }_e + \nu _\mu \,. \end{eqnarray} \tag{ 45 }$$

For the energy spectrum from weak decays, we follow Lipari et al. (2007, their Section IV). In this case, the decays are described by Equation (2). The functions dn^IT_a→b/dE_b simplify, in a frame where the parent a is ultra-relativistic, to dn^IT_a→b/dE_b = 1/E_a F_a→b(E_b/E_a). The functions F_a→b include the measured branching ratios in the possible final states and are given in Lipari et al. (2007, their Section IV). Note that π⁺ and π⁻ are initially produced in different ratios and produce muons with different helicities, described by the scaling functions with r_π = (m_μ/m_π)²:

$$\begin{eqnarray} F_{\pi ^+\rightarrow \mu _R^+}(x)=&F_{\pi ^-\rightarrow \mu _L^-}(x)=\frac{r_\pi (1-x)}{(1-r_\pi)^2x}\Theta (x-r_\pi)\qquad \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} F_{\pi ^+\rightarrow \mu _L^+}(x)=&F_{\pi ^-\rightarrow \mu _R^-}(x)=\frac{x-r_\pi }{(1-r_\pi)^2x}\Theta (x-r_\pi).\qquad \end{eqnarray} \tag{ 47 }$$

The muons decay further in a helicity-dependent way:

$$\begin{eqnarray} F_{\mu ^+\rightarrow \bar{\nu }_\mu }(x,h)&= &F_{\mu ^-\rightarrow \nu _\mu }(x,-h)\nonumber\\ &=&\left(\frac{5}{3}-3x^2+\frac{4x^3}{3}\right)+h\left(-\frac{1}{3}+3x^2-\frac{8x^3}{3}\right)\nonumber\\ \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} F_{\mu ^+\rightarrow \nu _e}(x,h)&= &F_{\mu ^-\rightarrow \bar{\nu }_e}(x,-h)\nonumber\\ &=&(2-6 x^2+4 x^3)+h(2-12 x+18x^2-8 x^3),\nonumber\\ \end{eqnarray} \tag{ 49 }$$

with h = 1 for right-handed and h = −1 for left-handed muons. It is therefore mandatory to distinguish four muon states μ⁺_L, μ⁺_R, μ⁻_L, and μ⁻_R as final states in order to account for the impact of muon polarization. The decay rate Γ^IT_a→b = Γ in Equation (2) (there is only one interaction type, which is decay) is just the inverse lifetime Γ = τ⁻¹ = τ⁻¹₀ m_a E⁻¹_a, where τ₀ is the rest-frame lifetime.

For kaons, the leading decay mode into muon and neutrino is treated in the same way as in Lipari et al. (2007) for the pion decays, i.e., with m_π → m_K. The branching ratio for this channel is about 63.5%. The second-most-important decay mode is K^± → π^± + π⁰ (20.7%). The other decay modes account for 16%, no more than about 5% each. Since interesting effects can only be expected in the energy range with the most energetic neutrinos, we only use the direct decays from the leading mode.

For protons accelerated in the jet, neutrons are produced by photohadronic interactions as described in Appendix B. Assuming that the neutrons escape from the acceleration region before they interact again, an additional neutrino flux from neutron decays is obtained. In this section, we show this neutrino flux separately. The beta decay describes the decay of the neutron into a proton, an electron, and an electron antineutrino. In the ultra-relativistic case, the mean fraction of the neutron energy going into the neutrino is χ ≈ 5.1 × 10⁻⁴; see Lipari et al. (2007). The neutrino spectrum is therefore obtained from the following equation:

$$\begin{equation} Q_{\bar{\nu }_e}(E_\nu)=\frac{1}{\chi } Q_n\left(\frac{E_\nu }{\chi }\right), \end{equation} \tag{ 50 }$$

with Q_n calculated from Equation (B6) in Appendix B.

We show the ν_μ neutrino spectra obtained from Sim-B and SOPHIA in Figure 10, where we also show the different contributions from different decay modes separately. Obviously, the SOPHIA and our combined spectra match very well, apart from the already discussed difference in the kinematics leading to some averaging in SOPHIA for large energies. Since the production of π⁺ dominates for initial protons, ν_μ in Figure 10 is most abundantly produced from pion decay; cf. Equation (44). However, the ν_μ from muon decays, coming from the π⁻ decay chain—cf. Equation (45)—is found at slightly higher energies and dominates for every high energy in the spectrum. For $\bar{\nu }_\mu$, the situation is exactly the opposite, but the final spectra look very similar. For very high energies, the SOPHIA spectrum is slightly higher than what one would expect because other decay modes (such as from neutral kaons) contribute, which we have not considered. Without synchrotron cooling, the contribution from kaon decays is small, however.

Zoom In Zoom Out Reset image size

In Figure 11, we show the ν_e (dashed) and $\bar{\nu }_e$ (solid) fluxes for our benchmarks. Obviously, the ν_e, coming from the π⁺ decay chain in Equation (44), dominate over the $\bar{\nu }_e$. However, if the neutrons produced by the photohadronic interactions escape from the source and then decay, they will lead to an additional neutrino flux shown by the thin gray curves (not included in the total $\bar{\nu }_e$ curves). Especially at very low energies, the $\bar{\nu }_e$ flux then dominates.

Zoom In Zoom Out Reset image size

5.3. Flavor and Neutrino–Antineutrino Ratios of the Neutrinos

In this section, we discuss the electron-to-muon neutrino flavor ratio (the ratio between the electron and muon neutrino fluxes) and the neutrino–antineutrino ratios at the source. The flavor composition at the source is primarily characterized by the flux ratio $(\bar{\nu }_e+ \nu _e)/(\bar{\nu }_\mu + \nu _\mu)$, since almost no ν_τ (or $\bar{\nu }_\tau$) are expected to be produced at the source. Because neutrino telescopes cannot distinguish neutrinos from antineutrinos in muon tracks or showers, this ratio is representative for the detection as well. From Equations (44) and (45), we can determine that this flavor ratio should be about 1/2 without kinematical effects, which we call the "standard assumption." On Earth, the three neutrino flavors then almost equally mix (in the ratio 1:1:1) through neutrino flavor mixing (Learned & Pakvasa 1995).

In Figure 12 we show the flavor ratio at the source for the GRB (left panel), AGN (middle panel), and BB (right panel) benchmarks. The standard assumption, 1/2, is marked by the horizontal lines. Our curves from Sim-B (thick solid neutrinos from pion and muon decays only) bend upward from this standard assumption, whereas the SOPHIA curves bend downward for large energies. This difference can be explained by a different implementation of the weak decays: if we do not take into account the spin state of the final muon (dotted curves h = 0), we can reproduce the SOPHIA results almost exactly with Sim-B. In fact, the effect of the helicity is larger than the details of the interaction model. The dashed curves include the effect of neutron decays (low energies) and kaon decays (high energies) into Sim-B. Especially for low energies, where $\bar{\nu }_e$ are abundantly produced by neutron decays, the curves deviate from the thick ones. This effect is strongest for the AGN benchmark, for which the standard assumption, 1/2, only approximately holds in a relatively narrow energy window. Note that the dashed curve in the left panel and all the other results for the GRB benchmark exactly match Lipari et al. (2007), where the weak decays are discussed in detail.

Zoom In Zoom Out Reset image size

If the Glashow resonance process $\bar{\nu }_e + e^- \rightarrow W^- \rightarrow {\rm anything}$ at around 6.3 PeV can be observed in a neutrino telescope, the neutrino–antineutrino ratios at the source may be relevant as well (all flavors at the source contribute to the $\bar{\nu }_e$ flux at the Earth through flavor mixing; Learned & Pakvasa 1995; Anchordoqui et al. 2005; Bhattacharjee & Gupta 2005). The neutrino–antineutrino ratio may be relevant to distinguish pγ interactions at the source, for which mostly π⁺ are produced, from pp interactions at the source, for which π⁺ and π⁻ are produced in almost equal ratios. Therefore, in Figure 13 we show the neutrino–antineutrino ratios at the source for electron neutrinos (left panels) and muon neutrinos (right panels). The electron neutrino–antineutrino ratios in the left panels depend on the ratio of π⁺ and π⁻ produced; see Equations (44) and (45). Our result matches SOPHIA very well, especially in the important energy range from the peak of the spectrum two decades down, apart from the discrepancy for high energies for the BB benchmark (upper right panel) which we discussed already as it is coming from the π⁺/π⁻ ratio. The deviation between SOPHIA and Sim-B can be up to 30%. After flavor mixing, the correction to the electron neutrino–antineutrino ratio at the Earth is at the level of 10%, much smaller than the effect on the flavor mix expected from pp interactions. The muon neutrino–antineutrino ratios in the right panels do not depend on the ratio of π⁺ and π⁻ produced, as in every pion decay the same number of muon neutrinos and antineutrinos is produced; see Equations (44) and (45). In this case, our results match SOPHIA and the standard prediction very well, and the effects of neutron and kaon decays are small in the absence of cooling.

Zoom In Zoom Out Reset image size

We have discussed simplified models for photo-meson production in cosmic accelerators. The main purpose of this simplification has been the definition of a photohadronic interaction model useful for efficient modern time-dependent AGN and GRB simulations, and for large-scale parameter studies, such as of neutrino flavor ratios. The major requirements have been listed in the introduction. For example, the secondaries (pions, kaons) are not to be integrated out, since their synchrotron cooling affects the neutrino flavor ratios.

We have first re-phrased the problem in terms of a two-dimensional "response function" to be folded with arbitrary photon and proton input spectra in order to compute the secondary fluxes. The key idea for our simplified models has been the factorization of this two-dimensional response function, which has allowed us to eliminate one of the integrations. In order to include kinematics that are as good as possible, we have then defined a discrete number of different interaction types with different characteristics based on the underlying physics of SOPHIA. The kinematics of more complicated interactions, such as direct production or multi-pion production, has been simulated by the introduction of multiple interaction types for each production channel. In a step-by-step fashion, we have simplified the resonance treatment in order to arrive at our simplified model Sim-B. It allows for the computation of pion spectra with only one integral, summed over about 10 interaction types, and can be easily adopted from our description. The extendibility of this approach has been demonstrated by showing how K⁺ fluxes can be added, once a suitable parameterization is found.

Similarly, the response function can be easily changed if new data are provided, new processes can be included, or systematics on the particle physics can be added. Of course, some effort has to be spent to find a suitable parameterization for each process.

We have demonstrated that our results match the output of SOPHIA sufficiently well. However, there are some differences due to the more refined kinematics treatment of SOPHIA, which effectively corresponds to one additional integration. For example, for very narrow spectral features, such as rapid cutoffs, the spectra are naturally more smeared out by SOPHIA. This is especially the case for high energetic interactions where multi-pion production is dominant, as can be seen in the BB (blackbody) benchmark. However, our approach is much simpler in the sense that the interaction rate and the folding with the proton spectra is automatically taken into account (cf. Appendix D). In addition, we have included the spin state of the final muon in the pion decays, as described in Lipari et al. (2007), which is not included in SOPHIA, which leads to differences in the neutrino flavor ratios; in fact, the electron-to-muon neutrino flavor ratio at the source is typically larger than 0.5, instead of smaller, as predicted without taking into account the spin state. In particular, we obtain ν_e:ν_μ:ν_τ≃ 1:1.85:0 for the GRB benchmark, 1:1.96:0 for the AGN benchmark, and up to 1:1.82:0 for the BB benchmark close to the spectral peaks. This means that especially the AGN benchmark behaves as pion beam in spite of the helicity dependence of the muon decays, whereas the BB benchmark shows the strongest deviation.

Since our approach has allowed us to discuss the leading interaction types separately, we have also shown differences from the Δ(1232)-resonance approximation (see also Mücke et al. 1999a, 2000b). For example, we have shown how multi-pion production modifies the characteristic shape of the GRB neutrino spectrum expected from the resonance approximation (which is flat in E² times the flux in the intermediate energy range). In fact, all of the resonances combined do not dominate in any significant part of the charged pion spectrum for our GRB, AGN, and BB benchmarks. In addition, the Δ(1232)-resonance approximation has been rendered insufficient for the computation of the (electron) neutrino–antineutrino ratios at the source because, by definition, no π⁻ is produced. We have also demonstrated from our general response function in model Sim-B, that for any input proton or photon spectrum the π⁺/π⁰ ratio at the source is significantly larger than one, as opposed to 1/2 from the Δ(1232)-resonance approximation. This implies that any neutrino flux based on the Δ(1232)-resonance approximation and normalized to photon flux observations using the charged-to-neutral pion ratio is underestimated by at least a factor of 2.4, where this factor is independent of the input spectra.

We conclude that our simplified model Sim-B allows for an efficient computation of pion and neutrino fluxes, including the necessary features for neutrino flux ratio discussion and the necessary efficiency for time-dependent simulations. It is sufficient for many purposes especially for power-law photon fields, but, of course, it cannot replace a full Monte Carlo simulation including full kinematics if high precision fits of existing data are required. In particular, it may turn out to be useful for time-dependent simulations and extensive parameter space studies using power-law spectra, including spectral breaks. However, in other cases, such as if the high-energy interactions with photon spectra with sharp peaks are very important, or there are anisotropic photon distributions, the Monte Carlo method may be the better choice.

We thank M. Maltoni, K. Mannheim, D. Meloni, and A. Reimer for useful discussions. S.H. acknowledges support from the Studienstiftung des deutschen Volkes (German National Academic Foundation). F.S. acknowledges support from Deutsche Forschungsgemeinschaft through grant SP 1124/1. W.W. acknowledges support from the Emmy Noether program of Deutsche Forschungsgemeinschaft, contract WI 2639/2-1. This work has also been supported by the Research Training Group GRK1147 "Theoretical astrophysics and particle physics" of Deutsche Forschungsgemeinschaft.

Here we follow the approach of Rachen (1996). As mentioned in Section 3.2, the direct production process is strongly backward peaked with respect to the produced pion. One possible approximation for the determination of χ is to assume that the average scattering angle is 180°. The energy fraction χ going into the pion as a function of _r is shown in Figure 14 as the curve cos θ = −1 (cf. Equation (16)). For comparison, the curves cos θ = 0 and cos θ = 1 are shown as well. For a more accurate representation of kinematics, we take the probability distribution of the Mandelstam variable t, as given in Mücke et al. (2000a), into account. For small |t|, it is given by

$$\begin{equation} \frac{dP}{dt}=b\frac{e^{b\,t(\cos {\theta _{\pi }})}}{e^{b\,t(-1)}-e^{b\,t(1)}} \end{equation} \tag{ A1 }$$

with b ≈ 12 GeV^-2. For example, for the interaction type T₁, the Mandelstam variable t is given by

$$\begin{eqnarray} t(\cos {\theta _{\pi }})&=& m_p^2-2 \frac{s+m_p^2}{2\sqrt{s}}\frac{s+m_p^2-m_\pi ^2}{2\sqrt{s}}\nonumber\\ &&-2\frac{s-m_p^2}{2\sqrt{s}}\sqrt{\left(\frac{s+m_p^2-m_\pi ^2}{2\sqrt{s}}\right)^2-m_p^2} \, \cos {\theta _{\pi }}\,.\nonumber\\ \end{eqnarray} \tag{ A2 }$$

Combining Equations (A1) and (A2), we obtain the average scattering angle as a function of the center-of-mass energy. Inserting the result into Equation (16), we obtain the average energy fraction χ going into the pion for direct one pion production. It is shown in Figure 14 as curve "〈cos θ〉" as a function of _r. Analogously, we compute the mean χ for direct two pion production by combining Equation (A1) with the variable t for the considered process.

Zoom In Zoom Out Reset image size

In the simplified model for direct production (cf. Section 4.3) we use the factorized response function (see Equation (27)) as for the simplified models of the other processes. Therefore, we have to choose an energy-independent, constant χ. Since the range of χ-values for direct production is wide, we divide the _r-range into three sections (for each of the interaction types T₁ and T_2a), such that it reproduces the results of the non-simplified model for typical power-law spectra for photons and protons in astrophysical sources, such as GRBs and AGNs. This approach corresponds to the thick gray curve "SIM" in Figure 14, where only the lower two interaction types are visible in the plotted energy range. Obviously, it is a good step function approximation to the curve "〈cos θ〉." The different sections of this step function correspond to our interaction types "L," "M," and "H."

A related issue to secondary particle production is the cooling or escape of the primaries due to the interaction process. We do not focus on the cooling or escape timescales in this study, but, for the sake of completeness, we demonstrate how they can be computed from the quantities presented for our simplified models in Section 4. If the primaries lose energy in an interaction, such as protons or neutrons in pion photoproduction p + γ → p + π⁰, this process can be interpreted as cooling, whereas if the primaries disappear, such as protons in p + γ → n + π⁺, it can be interpreted as escape. In the latter case, neutrons are re-injected into the system, which we will discuss below. The cooling rate t⁻¹_cool(E_p) = −E⁻¹_pdE_p/dt or escape rate t⁻¹_esc = −N⁻¹_pdN_p/dt for the species p (proton or neutron) due to the photohadronic interactions is, for constant K^IT, given by

$$\begin{eqnarray} t_{\mathrm{cool}}^{-1}(E_p) &=& \sum _{\mathrm{IT}} \Gamma _{p \rightarrow p}^{\mathrm{IT}}(E_p) \, K^{\mathrm{IT}} \,, \quad\nonumber\\ t_{\mathrm{esc}}^{-1}(E_p) &=& \sum _{\mathrm{IT}, \, p^{\prime } \ne p} \Gamma ^{\mathrm{IT}}_{p \rightarrow p^{\prime }}(E_p) \end{eqnarray} \tag{ B1 }$$

in terms of the quantities in Equation (2). Here, K^IT · E_p is the loss of energy per interaction; therefore, K^IT is often called "inelasticity." Note that if it is a function of the kinematical variables, such as the center-of-mass energy, it has to be folded into the calculation of the interaction rate in t⁻¹_cool. However, in Section 4, we have constructed the interaction types such that K^IT is a constant. For photo-pion production, the inelasticity can be related to the χ^IT_a→b in Equation (6) by

$$\begin{equation} K^{\mathrm{IT}} = \sum _{b \ne p} \chi _{p \rightarrow b}^{\mathrm{IT}} \, M_b^{\mathrm{IT}} \,, \end{equation} \tag{ B2 }$$

i.e., the energy loss of the nucleon equals the energy deposited in all interaction products (other than the initial nucleon). Note that the classification as cooling or escape also depends on if protons and neutrons are distinguished in the final state. In this appendix, we distinguish protons and neutrons. In addition, note that in astrophysical objects there may be other sources of cooling and escape to be taken into account, such as synchrotron cooling or escape from the production region.

The quantity needed for the computation of Equation (B1) is the interaction rate in Equation (4). Comparing Equation (4) with Equation (27), we find for our simplified models that the interaction rate for interaction type IT of the initial nucleon p is

$$\begin{equation} \Gamma ^{\mathrm{IT}}(E_p) = \int _{\frac{\epsilon _{\mathrm{th}} m_p}{2 E_p}}^{\infty } d \varepsilon \, n_\gamma (\varepsilon) \, f^{\mathrm{IT}}\left(\frac{E_p \varepsilon }{m_p} \right) \,, \end{equation} \tag{ B3 }$$

i.e., conveniently parameterized in terms of our simplified response function. Then the cooling and escape rates in Equation (B1) can be written in terms of the initial M_p or different nucleon $M_{p^{\prime }}$ multiplicity ($M_p+M_{p^{\prime }}=1$):

$$\begin{eqnarray} t_{\mathrm{cool}}^{-1}(E_p) &=& \sum _{\mathrm{IT}} M_p^{\mathrm{IT}} \, \Gamma ^{\mathrm{IT}}(E_p) \, K^{\mathrm{IT}} \,, \quad \nonumber\\ t_{\mathrm{esc}}^{-1}(E_p) &=& \sum _{\mathrm{IT}, \, p^{\prime } \ne p} M_{p^{\prime }}^{\mathrm{IT}} \, \Gamma ^{\mathrm{IT}}(E_p) \,. \end{eqnarray} \tag{ B4 }$$

The nucleon multiplicities are for the resonances in model A given in Table 1 for the interaction types in Table 3, for the resonances in model B in Table 4, for multi-pion production $M_{p=p^{\prime }}=0.69$ and $M_{p \ne p^{\prime }}=0.31$, and for direct production in Table 1 for the interaction types in Table 5. In resonance model A the inelasticities can be obtained from the χ^IT in Table 3 according to Equation (B2), i.e., K^IT 1,2,3 = ∑_πχ^IT_p→π M^IT_π, where the number of pions produced M^IT_π = 1 for ITs 1, 2a, and 2b, and M^IT_π = 2 for IT 3 (for IT 2, the pions in 2a and 2b are summed over). For resonance model B, they are given in Table 4, and for multi-pion production K^Multi-π ≃ 0.6, and direct production they are listed in Table 5 (here IT T_2b is not counted separately).

The re-injection rate p → p' for initial nucleons p can be obtained analogously to Equation (27) and Equation (28) from

$$\begin{equation} R^{\mathrm{IT}}(x,y)= \delta (x - (1-K^{\mathrm{IT}})) \, M_{p^{\prime }}^{\mathrm{IT}} \, f^{\mathrm{IT}}(y) \end{equation} \tag{ B5 }$$

as

$$\begin{eqnarray} Q_{p^{\prime }}^{\mathrm{IT}} &=& N_p \left(\frac{E_{p^{\prime }}}{1-K^{\mathrm{IT}}} \right) \frac{m_p}{E_{p^{\prime }}} \int _{\epsilon _{\mathrm{th}}/2}^\infty dy \, n_\gamma \left(\frac{m_p \, y \, (1-K^{\mathrm{IT}})}{E_{p^{\prime }}} \right) \nonumber\\ &&\times M_{p^{\prime }}^{\mathrm{IT}} f^{\mathrm{IT}}(y). \end{eqnarray} \tag{ B6 }$$

Note that double counting of the same interaction has to be avoided. In particular, interaction type 2 must not be counted twice.⁷ In addition, note that for multi-pion production Equation (38) based on the total cross section should be used in all cooling, escape, and re-injection rates.

Our benchmarks are given in the SRF. The benchmark for GRBs is taken from Lipari et al. (2007). The photon spectrum, a broken power law, is given by

$$\begin{equation} n_\gamma (\epsilon)=\left\lbrace \begin{array}{@{}ll@{}} {N}_\gamma \left(\frac{\varepsilon }{\mathrm{GeV}}\right)^{-1} \frac{1}{\mathrm{GeV \, cm^3}}& 0.2 \,\mathrm{eV}\, \le \varepsilon \le 1 \, \mathrm{keV}\\ \ms \,{N}_\gamma \, 10^{-6} \left(\frac{\varepsilon }{\mathrm{GeV}}\right)^{-2} \frac{1}{\mathrm{GeV \, cm^3}}& 1 \,\mathrm{keV}\, \le \varepsilon \le 300 \, \mathrm{keV} \end{array}\right. \end{equation} \tag{ C1 }$$

and the proton spectrum is given by

$$\begin{eqnarray} n_p(E_p)&=& {N_p} \, \left(\frac{E_p}{\mathrm{GeV}}\right)^{-2} \exp {\left[-\left(\frac{E_p}{6.9\cdot 10^8\,\mathrm{GeV}}\right)^2\right]}\nonumber\\ &&\times \frac{1}{\mathrm{GeV \, cm^3}} \quad E_p \ge 1\, \mathrm{GeV}. \end{eqnarray} \tag{ C2 }$$

Note that there are dimensionless normalization constants N_γ and N_p. The resulting neutrino spectrum is characterized by a wide maximum in E²Q_ν(E). This benchmark is designed to fit an average of the total distribution of GRB. We limit ourselves to this average, since taking all extreme cases of GRB would go beyond the scope of this paper.

The benchmark for AGNs is adopted from Mücke & Protheroe (2001). The photon spectrum, a power law with a sharp cutoff, is given by

$$\begin{eqnarray} &&n_\gamma (\epsilon)=\left\lbrace \begin{array}{@{}l@{}} {N}_\gamma \, \left(\frac{\varepsilon }{\mathrm{GeV}}\right)^{-1.6}\frac{1}{\mathrm{GeV \, cm^3}}\quad \ 10^{-3}\,\mathrm{eV} \le \varepsilon \le 140 \, \mathrm{eV}\\ \ms {N}_\gamma (1.4 \times 10^{-7})^{0.2}\\ \quad \times \left(\frac{\varepsilon }{\mathrm{GeV}}\right)^{-1.8}\frac{1}{\mathrm{GeV \, cm^3}}\quad 140 \,\mathrm{eV} \le \varepsilon \le 3.6 \, \mathrm{keV} \end{array}\right.\nonumber\\ \end{eqnarray} \tag{ C3 }$$

and the proton spectrum is given by

$$\begin{eqnarray} n_p(E_p)&=&{N_p} \, \left(\frac{E_p}{\mathrm{GeV}}\right)^{-2} \exp {\left[-\left(\frac{E_p}{2.6\times 10^9 \, m_p}\right)^2\right]}\nonumber\\ &&\times \frac{1}{\mathrm{GeV \, cm^3}} \quad E_p \ge 1\,\mathrm{GeV} \end{eqnarray} \tag{ C4 }$$

with m_p = 0.938 GeV. This benchmark is well in the range of usual parameters of the HBLs (a subclass of AGNs), which are the most interesting objects for state-of-the-art Air Cerenkov Telescopes.

The third benchmark, the most challenging one, is a high energetic BB spectrum, adopted from Mücke et al. (1999b). The photon spectrum of temperature 10 eV is given by

$$\begin{eqnarray} n_\gamma (\epsilon)&=& {N}_\gamma \, 1.318\times 10^{31}\left(\frac{\varepsilon }{\mathrm{GeV}}\right)^{2}\frac{1}{\exp \left[\frac{\varepsilon }{\mathrm{GeV}}\times 10^8\right]-1}\nonumber\\ &&\times \frac{1}{\mathrm{GeV \, cm^3}} \end{eqnarray} \tag{ C5 }$$

and the proton spectrum with a sharp cutoff is

$$\begin{equation} n_p(E_p)=\left\lbrace \begin{array}{@{}ll@{}} {N}_p \, \left(\frac{E_p}{\mathrm{GeV}}\right)^{-2}\frac{1}{\mathrm{GeV \, cm^3}}& 10^{6}\,\mathrm{GeV} \le E_p \le 10^{12} \, \mathrm{GeV}\\ \, 0\, & \mathrm{else.}\end{array}\right. \end{equation} \tag{ C6 }$$

The BB temperature is designed to fit usual BLR photons. It is the most challenging benchmark because the proton spectrum is not smeared out for the highest energies by an exponential cutoff and the photohadronic interactions are dominated by high energetic multi-pion production due to the peaked photon spectrum considered. Even though the benchmark is adopted from Mücke et al. (1999b), it does not exactly represent real physics. Thermal spectra are usually produced outside the shock and are therefore beamed. Beaming is negligible only for the production of cosmogenic neutrinos with CMB photons (Stecker 1979).

In SOPHIA, an injected proton of energy E_p is assumed to interact for sure, and the secondary particle (of type b) distribution dn_p→b/dE_b(E_p, E_b) is computed for a specific proton energy or a range of proton energies. This means that

$$\begin{equation} \int \frac{d n_{p \rightarrow b}}{dE_b} (E_p,E_b) dE_b = N_b \end{equation} \tag{ D1 }$$

is the number of secondary particles produced. The particle spectrum can then be computed using

$$\begin{equation} Q_b(E_b) = \int dE_p N_p(E_p) \Gamma _{p\rightarrow b} (E_p) \, \frac{d n_{p \rightarrow b}}{dE_b} (E_p,E_b) \,. \end{equation} \tag{ D2 }$$

Note that this formula is very similar to Equation (2), but not split up into different interaction types.

As is obvious from Equation (D2), the interaction rate as a function of proton energy is needed as additional input. We use Equation (4) with the total cross section as depicted in Figure 1 and parameterized in Section 3.3 to compute the interaction rate (the cross section is already summed over all interaction types). As far as the units and normalization are concerned, it is useful to specify all energies in GeV and all cross sections in μbarn = 10⁻³⁰ cm² (note, however, that in SOPHIA, photon spectra are always given in eV). In this case, the interaction rate carries units of 3 × 10⁻²⁰ N_γ s^-1, where N_γ is the dimensionless normalization of the photon spectrum in Appendix C. From Equation (D2), it is obvious that Q_b comes in units of 3 × 10⁻²⁰ N_γ N_p s^-1 GeV^-1 cm^-3 if N_p is given in units of GeV^-1 cm^-3 and carries the (dimensionless) normalization factor N_p. The same units apply to the results from our simplified models.

SOPHIA uses logarithmic energy spacing in E_p and E_b and provides the output on a discretized energy grid in the form dn_p→b/dx_b(E_p, x_b) with x_b = log(E_b/E_p), equally spaced in Δx_p with x_p = log(E_p/GeV) and in Δx_b. For the easiest data extraction, it is advisable to use Δx_p = Δx_b. Then, it is useful to rewrite Equation (D2) as

$$\begin{equation} Q_b\big(x_b^j\big) = \Delta x_p \sum _{i} \, N_p\big(x_p^i\big) \, \Gamma _{p\rightarrow b} \big(x_p^i\big) \, \frac{d n_{p \rightarrow b}}{dx_b} \big(x_p^i,x_b^j\big) \, 10^{-x_b^j} \,, \end{equation} \tag{ D3 }$$

where the last factor comes from switching to the logarithmic scale. With Equation (D3), the SOPHIA output can be used directly, collecting all entries for a specific x^j_b from all proton energy bins. Note, however, that the particular output format still has to be taken into account (SOPHIA first lists the bin range filled with data, then the filled bins, as a function of the proton energy).

For the test runs of the AGN and GRB spectra, we use a proton energy grid between 1 GeV and 10¹⁰ GeV with 100 bins, i.e., Δx_p = 0.1. In addition, we use 100 output bins with a step size Δx_b = 0.1. For the GRB benchmark, we use 100,000 trials per proton bin, for the AGN benchmark 25,000 trials. For the BB test run, we use a proton energy grid between 10⁶ GeV and 10¹² GeV with 60 bins, i.e., Δx_p = 0.1. In addition, we use 75 output bins with a step size Δx_b = 0.1 and 10,000 trials per proton bin.