# Simplified models for photohadronic interactions in cosmic accelerators

## Abstract

We discuss simplified models for photo-meson production in cosmic accelerators, such as Active Galactic Nuclei and Gamma-Ray Bursts. Our self-consistent models are directly based on the underlying physics used in the SOPHIA software, and can be easily adapted if new data are included. They allow for the efficient computation of neutrino and photon spectra (from decays), as a major requirement of modern time-dependent simulations of the astrophysical sources and parameter studies. In addition, the secondaries (pions and muons) are explicitely generated, a necessity if cooling processes are to be included. For the neutrino production, we include the helicity dependence of the muon decays which in fact leads to larger corrections than the details of the interaction model. The separate computation of the , , and fluxes allows, for instance, for flavor ratio predictions of the neutrinos at the source, which are a requirement of many tests of neutrino properties using astrophysical sources. We confirm that for charged pion generation, the often used production by the -resonance is typically not the dominant process in Active Galactic Nuclei and Gamma-Ray Bursts, and we show, for arbitrary input spectra, that the number of neutrinos are underestimated by at least a factor of two if they are obtained from the neutral to charged pion ratio. We compare our results for several levels of simplification using isotropic synchrotron and thermal spectra, and we demonstrate that they are sufficiently close to the SOPHIA software.

astroparticle physics, neutrinos, elementary particles, galaxies: active

## 1 Introduction

Photohadronic interactions in high-energy astrophysical accelerators are, besides proton-proton interactions, the key ingredient of hadronic source models. 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 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 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 electron 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 disc, and the CMB may serve as 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 AGN have been described already in previous works (Mannheim, 1993; Mücke & Protheroe, 2001; Aharonian, 2002) , as well as leptonic models have been discussed by Maraschi et al. (1992); 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); Rüger et al. (2010). A necessary prerequisite for a time-dependent hadronic modeling is the efficient computation of the photohadronic interactions. The on-line calculation via Monte Carlo simulations is not feasible, therefore a parametric description is the most viable way; 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 arising from the decays of both primary pions and secondary muons. However, it was pointed out in LABEL:~(Rachen & Mészáros, 1998; 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 . 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 LABEL:~(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 be also used for source identification, see, e.g., Xing & Zhou (2006); Choubey & Rodejohann (2009). Except from flavor identification, the differentiation between neutrinos and antineutrinos could be useful for the discrimination between and 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 at around  (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 accurate enough 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. (2000), which is partially based on Rachen (1996), or are in very simplified approaches computed with the -resonance approximation

 p+γ→Δ+→{n+π+1/3of all casesp+π02/3of all cases. (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 (-channel) production of pions. Kaon production is included as well by the simulation of the corresponding QCD processes (fragmentation of color strings). The treatment in Eq. (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 the efficient use in every step of modern time-dependent source simulations of AGNs and GRBs. The treatment in Eq. (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 to obtain 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 include the intermediate particles explicitely. 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:

• The model should predict the , , and fluxes separately, which is needed for the prediction of photon, neutrino, and antineutrino fluxes, and their ratios.

• 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.

• The particle physics properties should be transparent, easily adjustable and extendable.

• The model should be applicable to arbitrary proton and photon input spectra.

• 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).

• 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.

• 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 Sec. 2, we review the basic principles of the particle interactions, in particular, the necessary information (response function) to compute the meson photo-production for arbitrary photon and proton input spectra. In Sec. 3, we summarize the meson photo-production as implemented in the SOPHIA software, but we simplify the kinematics treatment. In Sec. 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 Sec. 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 App. B how the cooling and escape timescales, and how the neutron/proton re-injection can be computed in our simplified models. This study is supposed to be written for a broad target audience, including particle and astrophysicists.

## 2 Basic principles of the particle interactions

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 (per energy interval) of daughter particles of species and energy from the decay of the parent particle can be written as

 Qb(Eb)=∑IT∫dEaNa(Ea)ΓITa→b(Ea)dnITa→bdEb(Ea,Eb). (2)

Here is the differential spectrum of parent particles1 (particles per energy interval), and is the interaction or decay rate (probability per unit time and particle) for the process as a function of energy (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 describes the distribution (as a function of parent and daughter energy) of daughter particles of type per final state energy interval . 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, we are not interested in, are typically integrated out. If more than one daughter particle of the same species 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 and are typically given per volume, but here this choice is arbitrary since it enters on both sides of Eq. (2).

Sometimes it is useful to consider the interaction or decay chain without being interested in , for instance, for the decay chain . In this case, one can integrate out . This approach has, for instance, been used in Kelner & Aharonian (2008) for obtaining the neutrino spectra. Note, however, that the parent particles or may lose energy before they interact or decay, such as by synchrotron radiation. We do not treat such energy losses in this study explicitely, 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 Eq. (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 is identified with the proton or neutron, which we further on abbreviate with or (for proton or neutron). In this case, the interaction rate depends on the interaction partner, the photon, as

 ΓITpγ→p′b(Ep)=∫dε+1∫−1dcosθpγ2(1−cosθpγ)nγ(ε,cosθpγ)σIT(ϵr). (3)

Here is the photon density as a function of photon energy and the angle between the photon and proton momenta , is the photon production cross section, and is the photon energy in the nucleon/parent rest frame (PRF) in the limit . The interaction itself, and therefore and , is described in the SRF. The daughter particles 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 of the photon distribution. This limits this specific model to scenarios where we have seed photons produced in the shock rest frame (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 can be replaced by one over , we have

 ΓITpγ→p′b(Ep)=12m2pE2p∞∫ϵthmp2Epdεnγ(ε)ε22Epε/mp∫ϵthdϵrϵrσIT(ϵr). (4)

Here is the threshold below which the cross sections are zero. Note that corresponds to the available center of mass energy of the interaction as

 s(ϵr)=m2p+2mpϵr. (5)

In Eq. (4), the integral over takes into account that the proton and photon may hit each other in different directions. In the most optimistic case, (head-on hit), whereas in the most pessimistic case, (photon and proton travel in 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 uni-directional photons this would require one more integration.

In general, the function in Eq. (2) is a non-trivial function of two variables and . For photo-meson production in the SRF, the energy of the target photons is typically much smaller than the incident nucleon energy, and . In this case,2

 dnITa→bdEb(Ea,Eb)≃δ(Eb−χITa→bEa)⋅MITb. (6)

The function , which depends on the kinematics of the process, describes which (mean) fraction of the parent energy is deposited in the daughter species. The function 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 (if it changes at a certain threshold, one can define different interaction types below and above the threshold). For the relationship between and the inelasticity (fractional energy loss of the initial nucleon), see App. B. As we will discuss later, Eq. (6) is a over-simplified for more realistic kinematics, such as in direct or multi-pion production, since is, in general, a more complicated function of and the initial state properties. In this case, the distribution is not sufficiently peaked around its mean, and the -distribution in Eq. (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, we will in these case define different interaction types with different values of , simulating the broad energy distribution after the integration over the input spectra.

As the next step, we insert Eqs. (6) and (4), valid for photoproduction of pions, into Eq. (2), in order to obtain:

 Qb(Eb)=∞∫EbdEpEpNp(Ep)∞∫ϵthmp2Epdεnγ(ε)Rb(x,y) (7)

with

 x≡EbEp,y≡Epεmp, (8)

and the “response function”

 Rb(x,y)≡∑ITRIT(x,y)≡∑IT12y22y∫ϵthdϵrϵrσIT(ϵr)MITb(ϵr)δ(x−χIT(ϵr)). (9)

Here is the secondary energy as a fraction of the incident proton or neutron energy, corresponds to the maximal available center of mass energy3, and the is the fraction of proton energy deposited in the secondary. Note that and are typically functions of , if they only depend on the center of mass energy of the interaction. We will define our interaction types such that is independent of .

If the response function in Eq. (9) is known, the secondary spectra can be calculated from Eq. (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 , summed over various interaction types. We show in

the cross section as a function of 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 . 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 in Eq. (9) maps the same region in in different regions of , depending on the interaction type. This means that while a particular cross section may dominate for certain energies , 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.

From the particle physics point of view, 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.

## 3 Review of the photohadronic interaction model

The description of photohadronic interactions is based on Rachen (1996); Mücke et al. (2000), which means that the fundamental physics is 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 towards an analytical description for the response function in Eq. (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:

-resonance region

The dominant resonance process is the -resonance (at ; cf., Eq. (5) for the --conversion):

 p+γ\lx@stackrelΔ(1232)→p′+π. (10)

This process produces neutral (for ) or charged (for ) pions. For instance, for protons in the initial state, see Eq. (1).

Higher resonances

The most important higher resonance contribution is the decay chain

 γ+p\lx@stackrelΔ,N→Δ′+π,Δ′→p′+π′ (11)

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

 γ+p\lx@stackrelΔ,N→ρ+p′,ρ→π+π′. (12)
Direct production

The same interactions as in Eq. (10) or Eq. (11) (with the same initial and final states) can also take place in the -channel, meaning that the initial and nucleon exchange a pion instead of creating a (virtual) baryon resonance in the -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 only takes place in the -channel, whereas in the process the -channel is possible as well, because the photon can only couple to charged pions. The -channel reactions typically have a pronounced peak in the -distribution, whereas they are almost structureless in the momentum transfer distribution. On the other hand, -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. The process is in SOPHIA Mücke et al. (2000) 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: production (for protons in the initial state) with the decay channel leading to highest energy neutrinos.4 Note that at even higher energies, other processes, such as charmed meson production, may contribute as well.

We show in

the total photo-meson cross section as a function of , analogously to Mücke et al. (2000) (barn = 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 Eq. (9) for each interaction type, we need the kinematics, entering , the multiplicities , and the cross section .

### 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 with () in the SRF. The Lorentz factor between CMF and SRF is:

 γcm=Ep+Eγ√s≃Ep√s (13)

with the total CMF energy from Eq. (5). The energy of the pion in the SRF can then be written as:

 Eπ=γcmEcmπ(1+βcmπcosθcmπ)=EpEcmπ√s(1+βcmπcosθcmπ) (14)

with , and the pion energy, momentum and angle of emission in the CMF. Note that backward scattering of the pion means that . From Eq. (14) we can read off the fraction of energy of the initial going into the pion (which is equal to the inelasticity of for ). Since is a two body reaction, the energies of and pion are given by Particle Data Group (2008)

 Ecmπ=s−m2p+m2π2√s (15)

with given by Eq. (5). Using Eq. (15) in Eq. (14), we can calculate in Eq. (9) as

 χ(ϵr)=s(ϵr)−m2p+m2π2s(ϵr)(1+βcmπcosθπ). (16)

Indeed, is a function of if is a constant. As described in Rachen (1996), the average 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 . This approximation is only very crude. Therefore, we calculate the mean for direct production by averaging the probability distribution of the Mandelstam variable , as we discuss in App. 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 Eq. (16) for the kaon production (with the rough estimate ), where the replacements and have to be made. Similarly, one has for the first and second pion for the interaction in Eq. (11) (Rachen, 1996)

 χa(ϵr) =s(ϵr)−m2Δ+m2π2s(ϵr)(1+βcmπcosθπ), (17) χb(ϵr) =s(ϵr)−m2π+m2Δ2s(ϵr)m2Δ−m2p+m2π2m2Δ(1+βcmΔcosθΔ), (18)

where , and for the interaction Eq. (12)

 χ(ϵr) =12s(ϵr)−m2p+m2ρ2s(ϵr)(1+βcmρcosθρ)=K2 (19)

with . Note that in order to evaluate the -function in Eq. (9) (if one integrates over ), one also needs the derivative , 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 add up to the number of pions produced (one or two), the multiplicities for the nucleons in the final state to one. They are needed in Eq. (6) and Eq. (B1) (cooling and escape timescale in App. 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 ) or for the direct production (-channel). Note that different resonances will contribute with a similar pattern, such as through the interaction types , , , , etc.. In addition note that the branching ratios into certain resonances and from a resonance into the described channels are absorbed in the cross sections. Therefore, the total yield is always one (or two) pions per interaction in Table 1.

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. (2000). Most of the energy lost by the proton (, or ) 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

, 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 multiplicities times cross sections in Figs. 9 and 10 of Mücke et al. (2000) more accurately, we choose (), (), and () for initial protons (neutrons). Close to the threshold, the decreasing phase space for pion production requires a modification of the threshold for () from initial protons (neutrons). This corresponds to a vanishing multiplicity below the threshold, i.e., we assume that below only two pions are produced. We assume for the fraction of the proton energy going into one pion produced in the multi-pion channel . In addition, we choose and for initial protons (or and for initial neutrons) in accordance with Fig. 11 of Mücke et al. (2000) for high energies. As 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., Sec. 4.4.2).

### 3.3 Cross sections

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

.

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

 σITBW(ϵr)=s(s−m2p)24π(2J+1)BγBoutsΓ2(s−M2)2+sΓ2=BIToutsϵ2rσIT0(ΓIT)2s(s−(MIT)2)2+(ΓIT)2s (20)

with given by Eq. (5). Here , , and being the spin, the nominal mass, and the width of the resonance, respectively. We consider the energy-dependent branching ratios and the resonances shown in the Appendix B of Mücke et al. (2000); 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 to 4). For simplicity, we neglect the slight cross section differences between - and -interactions and use the values for protons, which implies that we do not take into account the N(1675) resonance.5 In Table 2 we list all considered resonances and their contributions to the different interaction types. To account for the phase-space reduction near threshold the function is introduced and multiplied on Eq. (20):

 (21)

with and taking the values shown in Table 2.

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

 σT1(ϵr) =92.7Pl(ϵr,0.152,0.25,2)+40exp(−(ϵr−0.29)20.002)−15exp(−(ϵr−0.37)20.002) (22) σT2(ϵr) =37.7Pl(ϵr,0.4,0.6,2) (23)

with

 (24)

where . The cross sections are given in barns, 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 formulas (cross sections in barns):

 σMulti−π1(ϵr) =80.3Qf(ϵr,0.5,0.1)s−0.34 (25) σMulti−π2(ϵr) =[1−exp(−ϵr−0.850.69)](29.3s−0.34+59.3s0.095)  if ϵr>0.85. (26)

We have compared our resonant and direct production multiplicities times cross sections and proton to neutron ratios with SOPHIA (cf., Figs. 9 to 11 in Mücke et al. (2000)), and they are in excellent agreement.

## 4 Simplified models

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

### 4.1 Factorized response function

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

 RIT(x,y)=δ(x−χIT)MITbfIT(y)withfIT(y)≡12y22y∫ϵthdϵrϵrσIT(ϵr). (27)

The part describes at what energy the secondary particle is found, whereas the part describes the production rate of the specific species as a function of . If only the number of produced secondary particles is important, it is often useful to show the effective .

As the next steps, we evaluate Eq. (7) with Eq. (27) by integrating over and by re-writing the -integral as -integral:

 QITb=Np(EbχIT)mpEb∞∫ϵth/2dynγ(mpyχITEb)MITbfIT(y). (28)

This (single) integral is relatively simple and fast to compute if the simplified response function together with 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 , and for suitable interaction types.

### 4.2 Resonances

Here we describe two alternatives to include 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.

#### Simplified model A (Sim-A)

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

 (29)

where and correspond to surface area (in barn-) and position (in GeV) of the resonance, and are process-dependent constants. Therefore, using Eq. (29), the -integral in Eq. (9) can be easily performed, and we obtain again the simplified Eq. (27) with

 fIT(y)=12y2ϵIT,0rBITout^ΓITΘ(2y−ϵIT,0r)andχIT≡χIT(ϵIT,0r). (30)

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

#### 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