Neutrinos from Cosmic Accelerators including Magnetic Field and Flavor Effects

Abstract

We review the particle physics ingredients affecting the normalization, shape, and flavor composition of astrophysical neutrinos fluxes, such as different production modes, magnetic field effects on the secondaries (muons, pions, and kaons), and flavor mixing, where we focus on 𝑝𝛾 interactions. We also discuss the interplay with neutrino propagation and detection, including the possibility to detect flavor and its application in particle physics, and the use of the Glashow resonance to discriminate 𝑝𝛾 from 𝑝𝑝 interactions in the source. We illustrate the implications on fluxes and flavor composition with two different models: (1) the target photon spectrum is dominated by synchrotron emission of coaccelerated electrons and (2) the target photon spectrum follows the observed photon spectrum of gamma-ray bursts. In the latter case, the multimessenger extrapolation from the gamma-ray fluence to the expected neutrino flux is highlighted.

1. Introduction

In addition to gamma-ray and cosmic ray instruments, neutrino telescopes, such as IceCube [1] or ANTARES [2], provide interesting data on the sources of the highest-energetic particles found in the universe, so-called “cosmic accelerators”; see [3–6] for reviews. In particular, neutrinos are a prominent way to search for the origin of the cosmic rays or to discriminate between leptonic and hadronic models describing the observed spectral energy distribution of photons. There are numerous possible sources, see [7] for an overview and [8] for the general theory. Very interesting extragalactic candidates for neutrino and cosmic ray production may be gamma-ray bursts (GRBs) [9] and active galactic nuclei (AGNs) [10–12]. The most stringent bounds for these sources, which are expected to be roughly uniformly distributed over the sky, so far come from IceCube, which has recently released data on time-integrated [13] and time-dependent [14] point source searches, GRB neutrino searches [15], and diffuse flux searches [16]. So far, no astrophysical neutrinos have been detected, which has been for a long time consistent with generic Waxman-Bahcall [17] and Mannheim-Protheroe-Rachen [18] bounds. However, data from IC40 and IC59, referring to the 40 and 59 string configuration of IceCube, respectively, start to significantly exceed these bounds, see [15, 16, 19], which is in tension with the corresponding neutrino production models, such as [9, 20, 21] for GRBs. For example, neutrino data may soon challenge the paradigm that GRB fireballs are the sources of the ultra-high-energy cosmic rays (UHECRs) [22]. For constraints to AGN models, see, for example, [23]. As a consequence, the age of truth has come for neutrino astrophysics, which is also the age of precision. Especially since data are available now, it is necessary to critically review the underlying assumptions from both the astrophysics and particle physics perspectives and to develop the models from rough analytical estimates into more accurate numerical predictions.

In this paper, we focus on the minimal set of particle physics ingredients for the neutrino production, which must be present in virtually all sources, using several specific examples. We focus on photohadronic (𝑝𝛾) interactions for the meson production, with the exception of Section 4.3. We do not only discuss the predicted neutrino flux, but also the flavor and neutrino-antineutrino composition at source and detector. The discussed effects include the following.(i) Additional pion production modes, such as 𝑡-channel (direct) and multipion production; see, for example, [24–30].(ii) Neutron and kaon production; see, for example, [29–36] (in general, there is also an additional contribution from charmed meson production, see [32, 33] for a detailed comparison).(iii) The cooling and decay of secondaries (pions, muons, and kaons); see, for example, [29, 30, 33, 34, 37–39].(iv) Flavor mixing and possible new physics effects; see [40] for a review.(v) The helicity-dependence of the muon decays; see, for example, [28, 38] (this effect has been discussed earlier in the context of atmospheric neutrinos, see, for instance, [41–43]).(vi) Spectral effects, such as the energy dependence of the mean free path of the protons, and their impact on the prediction; see, for example, [44, 45].(vii) The impact of the maximal proton energy on the neutrino spectrum; see, for example, [36].(viii) Deviations from the frequently used 𝐸−2𝜈 neutrino flux assumption; see, for example, [46].

While many of these effects have been studied elsewhere in the literature, we mainly show examples generated with the NeuCosmA (Neutrinos from Cosmic Accelerators) software in this paper to present them in a self-consistent way.

The structure of this paper is as follows: in Section 2, we give a simplified picture for the connection among neutrinos, cosmic rays, and gamma-rays. Then in Section 3, we review the minimal set of ingredients for neutrino production from the particle physics perspective. In Section 4, we discuss neutrino propagation and detection, including the possibility to detect flavor and the use of the Glashow resonance, where we illustrate how new physics can be tested in the neutrino propagation in Section 4.4. We furthermore present two specific applications: a generic AGN-like model in Section 5 and a model for GRBs in Section 6, where the main difference is the model for the target photons. Then we finally summarize in Section 7.

2. Neutrinos and the Multimessenger Connection

Here we outline a simplified picture of the neutrino or cosmic ray source, as often used in the literature, whereas we add extra ingredients in the next section. In this approach, charged mesons originate from 𝑝𝑝 or 𝑝𝛾 interactions, where we focus on 𝑝𝛾 (photohadronic) interactions in this work; see, for example, [47, 48] for 𝑝𝑝 interactions, which may be dominant for particular source classes, such as supernova remnants. In the simplest possible picture, charged pions are produced by the Δ(1232)-resonance𝑝+𝛾⟶Δ+⟶{{{{{{{𝑛+𝜋+13ofallcases,𝑝+𝜋023ofallcases.(2.1)

While this process is not sufficient for state-of-the-art models for neutrino production, it is very useful to illustrate a few qualitative points common to many cosmic ray and neutrino production models. The protons on the l.h.s. of (2.1) are typically assumed to be injected into the interaction volume with an (𝐸′𝑝)−𝛼 spectrum (here primed parameters refer to the shock rest frame (SRF), whereas unprimed parameters to the observer’s frame) coming from Fermi shock acceleration, where 𝛼∼2. They interact with the photons on the l.h.s. of (2.1) with energy 𝜀′∼(0.2−0.3)GeV/𝐸′𝑝. While the assumptions for the injected protons are similar for most models (except from the minimal and maximal energies), the target photons are typically described in a model- and source-dependent way, such as by(1) synchrotron emission from coaccelerated electrons or positrons,(2) thermal emission, such as from an accretion disk,(3) a more complicated combination of radiation processes,(4) an estimate inferred from the gamma-ray observation,

just to name a few examples. While any realistic simulation of a particular source will imply option (3), the other options typically rely on fewer parameters and may be good approximations in many cases. In particular, a reliable prediction for the photon density in the source may be obtained from the gamma-ray observation, option (4), if the photons can escape. In fact, we will use option (1) in Section 5 and option (4) in Section 6.

After an interaction between proton and photon, the particles on the r.h.s. of (2.1) are produced with the given branching ratios. The neutrinos then originate from 𝜋+ decays via the decay chain𝜋+⟶𝜇++𝜈𝜇,𝜇+⟶𝑒++𝜈𝑒+𝜈𝜇,(2.2) where in this standard picture 𝜈𝑒∶𝜈𝜇∶𝜈𝜏 are produced in the ratio 1 : 2 : 0 if the polarities (neutrinos and antineutrinos) are added. In addition, high-energy gamma-rays are produced by𝜋098.8%−−−−−→𝛾+𝛾.(2.3)

These are typically emitted from the source at lower energies due to electromagnetic cascades, in addition to gamma-rays escaping from the interaction volume (the ones contributing on the l.h.s. of (2.1)).

From (2.1), we can also illustrate the production of cosmic ray protons, ignoring for the moment that the composition of cosmic rays may be heavier at high energies [49]. First of all, some of the protons injected into the interaction volume on the l.h.s. of (2.1) may escape, leading to cosmic ray production. However, even if the protons are magnetically confined, the neutrons on the r.h.s. of (2.1), which are electrically neutral, can easily escape if the source is optically thin to neutron escape. After decay (typically outside the source)𝑛⟶𝑝+𝑒−+𝜈𝑒,(2.4) they lead to cosmic ray flux and an additional 𝜈𝑒 neutrino flux which is an unavoidable consequence of the interactions in (2.1). The cosmic ray protons with energies above 6⋅1019 eV interact with the cosmic microwave background (CMB) photons by (2.1), leading to the so-called Greisen-Zatsepin-Kuzmin (GZK) cutoff [50, 51]. However, according to (2.1), charged pions are produced in these interactions as well, which means that an additional neutrino flux should come with that, which is often called “cosmogenic neutrino flux.”

In summary, the photohadronic interaction in (2.1) offers a self-consistent picture for a cosmic ray source, with a possible connection among cosmic ray, neutrino, and gamma-ray escape. In specific models, however, that does not mean that a large neutrino flux is guaranteed for every cosmic ray source. For instance, the interaction rate for the process in (2.1), which depends on the photon density, may be low.

3. Simulation of Neutrino Sources

Here we give a more detailed generic picture of the simulation of neutrino sources from the particle physics perspective with the minimal set of ingredients, without using a specific model. A flowchart summarizing the contents of this section, which can be followed during the reading, is given in Figure 6.

3.1. Photohadronic Interactions

In order to describe the processes within an interaction volume (one zone in the simplest case), two kinds of spectra are needed: 𝑄′(𝐸′) (in units of GeV⁻¹ cm⁻³ s⁻¹) describes the number of particles injected or ejected per volume element and energy interval, and 𝑁′(𝐸′) (in units of GeV⁻¹ cm⁻³) describes the particle density per energy interval. The secondary meson injection rate 𝑄′𝑏(𝐸′𝑏) for a pion or kaon species 𝑏 produced in photohadronic interactions is given by (following [28])𝑄′𝑏(𝐸′𝑏)=∫∞𝐸′𝑏𝑑𝐸′𝑝𝐸′𝑝𝑁′𝑝(𝐸′𝑝)∫∞0𝑐𝑑𝜀′𝑁′𝛾(𝜀′)𝑅𝑏(𝑥,𝑦).(3.1)

Here 𝑥=𝐸′𝑏/𝐸′𝑝 is the fraction of energy going into the secondary, 𝑦≡(𝐸′𝑝𝜀′)/𝑚𝑝 (here 𝑦 can be related to the center-of-mass energy by 𝑠=𝑚2𝑝+2𝑚𝑝(1−cos𝜃𝑝𝛾)𝑦, where 𝜃𝑝𝛾 is the angle between proton and photon momentum; 𝜃𝑝𝛾=𝜋 corresponds to heads-on collisions), and 𝑅𝑏(𝑥,𝑦) is the “response function.” If many interaction types are considered, the response function can be quite complicated. However, if it is known from particle physics, (3.1) can be used to compute the secondary injection for arbitrary proton and photon spectra. The important point here is that the secondary production depends on the product normalization of the proton density 𝑁′𝑝(𝐸′𝑝) and the target photon density 𝑁′𝛾(𝜀′) within the interaction volume. Thus a higher proton density can be compensated by a lower photon density, and vice versa (strictly speaking, this degeneracy only holds in the absence of any other radiation process, e.g., inverse Compton scattering depends on the photon density individually, which means that the input spectral shapes (especially 𝑁′𝛾 in (3.1)) will be modified if that process contributes significantly). Another implication of (3.1) is that the secondary production depends on the densities within the source 𝑁′, not the injection rates 𝑄′. Of course, one cannot look into the source, but can only observe cosmic messengers escaping from the source. As we will demonstrate later, the observed/ejected photon or cosmic ray spectrum 𝑄′ is only directly representative for the corresponding density spectrum within the source 𝑁′ if “trivial” escape is the leading process, that is, 𝑄′=𝑁′/𝑡′esc with 𝑡′esc∼𝑅′/𝑐 and 𝑅′ the size of the interaction region. For this section, (3.1) is used as a starting point for the computation of the neutrino fluxes, where we do not discuss the origin of spectral shape and normalization of 𝑁′𝑝 and 𝑁′𝛾. In practice, typically an (𝐸′𝑝)−2 injection spectrum is assumed for the protons, as mentioned above, where the maximal energy is limited by synchrotron and adiabatic losses. The photon density may be a consequence of a complicated interplay of radiation processes. In either case, the derivation of these densities depends on the model, and we will show several examples in Sections 5 and 6.

Compared to a numerical approach, (2.1) has limitations to describe the meson production. First of all, additional pion production modes contribute, such as higher resonances, direct (𝑡-channel), and multipion production, which will also lead to 𝜋− production (cf. Figure 6). These are not as easy to describe as the Δ-resonance because of different shapes of the cross sections and more complicated kinematics. The Monte Carlo package SOPHIA [25] can deal with these interactions. In order to increase the efficiency, often parameterizations of SOPHIA are used, such as [28, 52]. In the following, we use [28] (model Sim-B), because the secondary muons and pions are needed explicitly. We show the impact of the resonant production (including higher resonances) on 𝜋+ production in Figure 1 for a typical GRB, an AGN, and a (10 eV) black body (BB) target photon field. As one can read off from this figure, the resonances always give a reasonable first estimate for the actual pion production, but quantitatively they only dominate at the breaks. In addition, multipion processes can change the spectral shape significantly, such as for the GRB example, which is a consequence of the cross section dependence on the center-of-mass energy. As a further limitation, note that (2.1) does not describe kaon production and subsequent decay into neutrinos, where the leading modes are given by𝑝+𝛾⟶𝐾++Λ⧵Σ,(3.2)𝐾+⟶𝜇++𝜈𝜇.(3.3)

(a)

(b)

(c)

Figure 1 

Contributions of different 𝜋+ production modes to the total spectrum for a typical GRB, an AGN, and a (10 eV) black body (BB) target photon density, figure taken from [28].

The branching ratio for the leading channel in (3.3) is about 64%. The second-most-important decay mode is 𝐾±→𝜋±+𝜋0 (20.7%). The other decay modes account for 16%, no more than about 5% each. Because 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.

In the literature, the Δ-resonance approximation in (2.1) is, even in analytical approaches, typically not taken literally. For example, a simple case is the approximation by Waxman and Bahcall [9] (“WB Δ-approx.”), for which one can write the response function as𝑅𝜋±(𝑥,𝑦)=0.5×𝛿(𝑥−0.2)×500𝜇barn×{{{{{{{{{{{{{{{{{02𝑦<0.2GeV,1−(0.2GeV2𝑦)20.2GeV≤2𝑦<0.4GeV,0.12GeV2(2𝑦)22𝑦≥0.4GeV,(3.4)which implies that charged pions are produced in 50% of all cases, and these take 20% of the proton energy. In addition, the width of the Δ-resonance is taken into account (note that this description is slightly more accurate than [9], which uses an additional integral approximation). In fact, this function peaks at 2𝑦≃0.4 GeV, which is higher than the threshold for photohadronic interactions—and it is even a little bit higher in the numerical calculation. The reason is that for the threshold often head-on collisions are assumed (𝜃𝑝𝛾=𝜋), whereas these only contribute a small part to the total number of interactions. Using (3.4) in (3.1) and rewriting the integral over 𝜀′ in one over 𝑦, it is easy to show that for power law spectra 𝑁′𝑝∝(𝐸′𝑝)−𝛼 and 𝑁′𝛾∝(𝜀′)−𝛽𝑄′𝜋(𝐸′𝜋)∝(𝐸′𝜋)−𝛼+𝛽−1.(3.5)

This means that the pion spectral index depends on both the proton and photon spectra, where it is inversely proportional to the photon spectral index. As we will see below, the neutrino spectrum follows the pion spectrum, which means that the assumption of an 𝐸−2𝜈 spectrum for the neutrinos, as it is often used in data analysis and many models in the literature, is only a valid assumption for 𝛽≃1—which is roughly observed for GRBs below the break. On the other hand, if the target photons come from synchrotron emission, such a hard photon spectrum is not possible, and 𝛽≳3/2 may be more plausible if the electrons are injected with a spectral index similar to the protons. As a consequence, the neutrino spectrum becomes harder. In addition, multipion processes in the photohadronic interactions will act in the same direction and make the neutrino spectrum even harder, compare Figure 1; see also [29] for a detailed comparison between the approximation in (3.4) and the numerics. Note that for 𝑝𝑝 interactions with “cold” (non-relativistic) protons, the 𝐸−2𝜈 assumption may be plausible [47]. We discuss the implications of the 𝐸−2𝜈 assumption for the detector response in Section 4.

3.2. Decays of Secondaries

The weak decays of pions and muons are described in detail in [38]. In general, in case of ultra-relativistic parents of type 𝑎, the distribution of the daughter particle of type 𝑏 takes a scaling form in order to obtain for the energy spectra𝑄′𝑏(𝐸′𝑏)=∑𝑎∫∞𝐸′𝑏𝑑𝐸′𝑎𝑁′𝑎(𝐸′𝑎)(𝑡′dec)−11𝐸′𝑎𝐹𝑎→𝑏(𝐸′𝑏𝐸′𝑎),(3.6) summed over all parent species. The functions 𝐹𝑎→𝑏 for pion-, kaon-, and helicity-dependent muon decays can be read off from [38] (Section 4). Consider the simplified case of a 𝛿-function for 𝐹𝑎→𝑏. For instance, for decays of neutrons in (2.4), one may approximate𝐹𝑛→𝜈𝑒=𝛿(𝐸′𝜈𝐸′𝑛−𝜒𝑛),(3.7) with 𝜒𝑛=5.1×10−4. If the neutrons do not interact, 𝑄′𝑛=𝑁′𝑛(𝑡′dec)−1 (see below), and we find for the neutrino injection  (3.6)𝑄′𝜈𝑒(𝐸′𝜈)=1𝜒𝑛𝑄′𝑛(𝐸′𝜈𝜒𝑛).(3.8)

In this case, the neutrino spectrum follows the neutron spectrum. If the neutrinos originate from pion decays, the neutrino injection follows the pion injection spectrum by similar arguments.

Since the decays of muons are helicity dependent, that is, 𝐹𝑎→𝑏 is different for left- and right-handed muons, it is necessary to keep track of these two species separately (cf. Figure 6). Although the effects of the helicity-dependent muon decays on the fluxes are probably small, the flavor composition is slightly affected, depending on the parameters of the source. We illustrate this effect in Figure 2 for the GRB, AGN, and black-body examples in Figure 1. In this figure, the horizontal line corresponds to the standard assumption, that is, neutrinos being produced in the flavor composition 𝜈𝑒∶𝜈𝜇∶𝜈𝜏 of 1 : 2 : 0. Including the scaling of the secondary decays, the dotted (blue) curves are obtained if the helicity of the muons is averaged over. From the comparison with the light (red) solid curves, it is clear that this assumption is implemented in SOPHIA. On the other hand, as pointed out in [38], keeping track of the muon helicity slightly changes the flavor composition, see dark (black) solid curves, which are significantly different from the standard assumption and the helicity-averaged version. However, it is also clear from Figure 2 that the deviation from the standard prediction depends on the input spectra, as it is smaller for the AGN than the GRB example. The dashed (green) curves show the contributions of the neutron and kaon decays, which affect the flavor composition at very low and high energies, respectively.

(a)

(b)

(c)

Figure 2 

Flux ratio at the source between electron and muon neutrinos as a function of the neutrino energy, for the same objects as in Figure 1. No secondary cooling is included in this figure, figure taken from [28].

3.3. Cooling of Secondaries

In order to describe the cooling of the secondary pions, muons, and kaons, we use the steady state approach; that is, we do not allow for an explicit time-dependence since the statistics of neutrino observations is typically expected to be low. The steady state equation for the particle spectrum, assuming continuous energy losses, is given by𝑄′(𝐸′)=𝜕𝜕𝐸′(𝑏′(𝐸′)𝑁′(𝐸′))+𝑁′(𝐸′)𝑡′esc,(3.9) with 𝑡′esc(𝐸′) the characteristic escape time, 𝑏′(𝐸′)=−𝐸′𝑡′−1loss with 𝑡′−1loss(𝐸′)=−1/𝐸′𝑑𝐸′/𝑑𝑡′ the rate characterizing energy losses. This differential equation balances the particle injection on the l.h.s. with energy losses and escape on the r.h.s. of the equation. Note again that the steady density 𝑁′(𝐸′) is needed for the photohadronic interactions in (3.1), not the actual injection spectrum. In addition, note that if there are no energy losses (𝑏′=0), one has immediately 𝑄′(𝐸′)=𝑁′(𝐸′)/𝑡′esc from (3.9), which we have already used above. If decay is the dominant escape mechanism, one finds 𝑄′(𝐸′)∝𝑁′(𝐸′)/𝐸′.

While the primary proton and photon spectra in (3.1) could be affected by a number of radiation processes, for the neutrino fluxes and flavor compositions, at least the processes of the secondaries (pions, muons, kaons) are important. We illustrate the synchrotron cooling and decay rates for pions, muons, and kaons as a function of energy in Figure 3. As one can easily see in the figure, for any species, decay dominates at low energies, while synchrotron cooling dominates at high energies. Other cooling or escape processes are often subdominant, as illustrated by the gray-shaded region for an adiabatic cooling components in GRBs. The two curves meet at a critical energy 𝐸′𝑐 for each species, which is different depending on the particle physics parameters. As a consequence of (3.9), the corresponding steady spectra 𝑁′ are loss-steepened by two powers above𝐸′𝑐=√9𝜋𝜖0𝑚5𝑐7𝜏0𝑒4𝐵′2,(3.10) where synchrotron and decay rates are equal. These critical energies depend on the particle physics properties of the parent, that is, the mass 𝑚 and the rest frame lifetime 𝜏0, and the magnetic field 𝐵′ as the only astrophysical parameter. It is therefore a very robust prediction and might allow for the only direct measurement of 𝐵′. Rescaling the magnetic field shifts the critical energies by a constant amount on the horizontal axis, but does not change the spacing in the logarithmic picture. In Figure 3, we also show the estimated range for adiabatic cooling as additional cooling component (shaded region), which may have some impact, especially on the muons, in extreme cases. In these cases, the height of the spectral peaks will be somewhat reduced, but the qualitative picture does not change.

Figure 3 

Synchrotron cooling and decay rates for pions, muons, and kaons as a function of energy for a magnetic field of 𝐵′≃300 kG. The gray-shaded region shows the estimated range for the escape (or adiabatic cooling) rate for GRBs for a variability timescale ranging from 0.01 s to 1 s, figure taken from [30].

We show in Figure 4 the consequences for the spectral shape at one GRB example. Here the injection spectrum of electron neutrinos Figure 4(a) and muon neutrinos Figure 4(b) is shown, including the individual contributions from the parents. First of all, one can read off from this figure that in the case of no losses (dashed curves) the spectral shapes of all contributions are very similar, and the neutrino fluxes add in a trivial manner. A change of the primary spectra 𝑁′𝑝 and 𝑁′𝛾 in (3.1) may change the shape of the dashed curves, but almost in the same way for all curves. If the synchrotron losses are switched on (solid curves), the spectral split predicted by (3.10) (see also Figure 3) among the neutrino spectra coming from different parent species can be clearly seen in Figure 4(b). One can also see a small pile-up effect coming from the muon decays; that is, a small region where the cooled muons coming from higher energies pile up and dominate, and lead to a higher flux than in the “no losses” case. In Figure 4(a), only two spectra are shown, since only muon or neutron decays may produce the electron flavor. Because of the very small 𝜒𝑛 in (3.8), the neutron decays only show up at low energies. Since the decay of pions, muons, or kaons, which dominate in different energy ranges, leads to different neutrino flavor compositions, the flavor composition at the source will be changed as a function of energy. In the literature, often the following source classes in terms of 𝜈𝑒∶𝜈𝜇∶𝜈𝜏 are distinguished (see [53] for a summary, including other options leading to a neutron beam or muon beam-like source).

(a)

(b)

Figure 4 

Injection spectra for electron neutrinos (a) and muon neutrinos (b) at the source for a GRB example, where the individual contributions from the parents are shown (neutrinos and antineutrinos added), figure taken from [30].

Pion Beam Source
Pion decays and muon decays equally contribute; flavor composition 1 : 2 : 0.

Muon Damped Source
Strong muon cooling, which means that pion decays dominate; flavor composition 0 : 1 : 0.

Muon Beam Source
Pile-up muons dominate; flavor composition 1 : 1 : 0.

Neutron Beam Source
Neutron decays dominate; flavor comp. 1 : 0 : 0.

Undefined Source
Several of these processes compete; flavor composition 𝑋∶1−𝑋∶0 (𝑋∉{0,1/3,1/2,1}).

Note that in none of these cases 𝜈𝜏 are produced in significant quantities, which means that it is sufficient to use the ratio between the 𝜈𝑒 and 𝜈𝜇 production at the source, sometimes called “flavor ratio.” From the preceding discussion, it is clear that the above classifications can only hold in specific energy ranges. However, all of these sources can be recovered in a numerical simulation. In Figure 5, four examples for energy-dependent flavor ratios at the source are shown for different parameter sets, where in this model 𝑁′𝛾 is generated by synchrotron losses of coaccelerated electrons. In this figure, also the relevant flavor ratio ranges for the different sources introduced above are shown. Figure 5(b) shows the classical pion beam source, which is typically found for low magnetic fields. Nevertheless, the contribution of neutron decays at low energies can be clearly seen. Figure 5(d) shows a pion beam evolving in a muon damped source at high energies, as, for example, in [37]. Figure 5(a) depicts a muon beam to muon damped source. In this case, the cooled muons pile up at lower energies, where the muon decays dominate. And Figure 5(c) shows an undefined source, where several processes compete.

(a)

(b)

(c)

(d)

Figure 5 

Energy-dependent flavor ratio at the source for several selected parameter sets (test points “TP”; cf. Figure 10) for a model where 𝑁′𝛾 is generated by synchrotron losses of coaccelerated electrons. The (dark) gray-shaded areas mark the regions where the fluxes peak, figure taken from [34].

Figure 6 

Flowchart describing the model (in SRF). The functions 𝑄′(𝐸) denote (injection) spectra per time frame [(GeV cm3 s)−1] and 𝑁′(𝐸) steady spectra [(GeV cm3)−1] derived from the balance between injection and losses or escape. Dashed arrows stand for solving the steady state differential equation (3.9), figure taken from [30].

Of course, not only the secondaries are affected by synchrotron (or adiabatic) losses, but also the primary protons. Depending on the model, one can use these losses to determine the maximal proton energy, or one can put in the maximal proton energy by hand. From Figure 4, it is interesting to discuss the impact of the maximal proton energy on the neutrino fluxes. In the “no losses” case (dashed curves), the maximal neutrino energy is directly determined by the maximal proton energy 𝐸′𝜈,max≃0.05𝐸′𝑝,max. There is, however, one exception: the neutrinos from neutron decays are limited by 𝐸′𝜈,max≃10−4𝐸′𝑝,max, cf., (3.8). If the synchrotron losses are switched on (solid curves), the neutrino spectrum from neutron decays still follows the proton spectrum, since the neutrons are electrically neutral, whereas the maximal neutrino energies for the other production modes are determined by (3.10). As a consequence, the neutron decay spectrum strongly depends on the assumptions for 𝐸′𝑝,max, whereas the other spectral shapes are entirely unaffected by 𝐸′𝑝,max as long as 𝐸′𝑝,max≳6𝐸′𝑐 for kaons (see (3.10), factor six from kaon production and decay kinematics). It is therefore not surprising that for strong enough magnetic fields one can find parameter sets for which the neutron decays dominate (cf. [35, 36]), but one should keep in mind that this depends on the assumptions for the maximal proton energy (and the inclusion of multipion processes, etc., which may mask this effect).

3.4. Transformation into Observer’s Frame

The transformation of the injection spectrum of the neutrinos 𝑄′𝜈𝛼 from the source to the observable flux 𝜙𝛽 of 𝜈𝛽 (in units of GeV⁻¹ cm⁻² s⁻¹) at the Earth is given by𝜙𝛽=∑𝛼=𝑒,𝜇,𝜏̂‌𝑁𝑃𝛼𝛽(1+𝑧)24𝜋𝑑2𝐿𝑄′𝜈𝛼,𝐸𝜈=Γ1+𝑧𝐸′𝜈,(3.11) where a simple Lorentz boost Γ is used (instead of a viewing angle-dependent Doppler factor). Here 𝑃𝛼𝛽 is the transition probability 𝜈𝛼→𝜈𝛽, discussed in Section 4, and ̂‌𝑁 is a (model-dependent) normalization factor. For example, if an isotropically emitting spherical zone is boosted with Γ towards the observer, then ̂‌𝑁=(4/3)𝑅′3𝜋Γ2 since the emission is boosted into a cone with opening angle 1/Γ. For a relativistically expanding fireball, it is simpler to perform the transformation in a different way, see Section 6. Furthermore, 𝑑𝐿(𝑧)=(1+𝑧)𝑑com(𝑧) is the luminosity distance, and 𝑑com(𝑧) is the comoving distance. From (3.11), one can read off that the redshift dependence of the neutrino luminosity scales as 𝐸2𝜈𝜙∝1/𝑑2𝐿 independent of the model (for fixed 𝑄′𝜈𝛼), which is expected. Note that in (3.11) the neutrino and antineutrino fluxes are often added if the detector cannot distinguish these.

3.5. Summary of Ingredients, and Limitations of the Approach

We summarize the generic neutrino production in Figure 6. As one can see in this flowchart, the starting points are the proton and photon densities within the source. Once these (and 𝐵′) are fixed, the rest is just a particle physics consequence. Therefore, for the computation of specific neutrino fluxes, the main effort is actually to determine 𝑁′𝑝(𝐸′𝑝), 𝑁′𝛾(𝜀′), and 𝐵′.

Of course, there are some processes not taken into account in this picture, which may add to the ingredients discussed above for specific source classes. For instance, secondary neutrons, produced in the photohadronic interactions, may interact again if the source is optically thick to neutrons, the secondary pions, muons, and kaons may be reaccelerated [54], synchrotron photons of the secondaries may add to 𝑁′𝛾, and so forth. In addition, the neutrino spectrum may be more complicated in multizone models, since these naturally allow for more freedom.

From a particle physics perspective, additional kaon and charmed meson production modes may be added, and also the secondaries may interact again (see [32, 33]). However, compared to analytical computations, the numerical approach described in this section already takes into account the secondary cooling in a self-consistent way, additional neutrino production modes can be easily included, and the full energy dependencies can be accounted for. For example, it has been demonstrated in [44] that all necessary ingredients to reproduce the analytical GRB fireball (neutrino) calculations in [9, 15, 20, 21] are contained. It should represent the minimal set of of ingredients for neutrino production which are present in every source in the spirit of constructing the simplest possible model first. Of course, if 𝐵′ is small, the secondary cooling effects will be small as well, which is automatically included.

4. Neutrino Propagation and Detection

In this section, we discuss several aspects of neutrino propagation and detection from the theoretical perspective.

4.1. Neutrino Propagation and Observables

It is well known that neutrinos may change flavor from the production to the detection point. While this phenomenon is in general described by neutrino oscillations, astrophysical neutrinos are typically assumed to suffer from decoherence over very long distances. This means that effectively (in most practical cases) only flavor mixing enters the astrophysical neutrino propagation (see [55] for a more detailed discussion). In that case, 𝑃𝛼𝛽 in (3.11) becomes𝑃𝛼𝛽=3∑𝑖=1|𝑈𝛼𝑖|2|𝑈𝛽𝑖|2,(4.1) for three active neutrinos, where 𝑈𝛼𝑖 are the usual PMNS mixing matrix elements in the standard parameterization; see, for example, [56] for recent values of the mixing angles. This implies that neutrino oscillations, that is, the Δ𝑚2𝐿/𝐸-dependence, are averaged out. An initial flavor composition 𝜈𝑒∶𝜈𝜇∶𝜈𝜏 of 1 : 2 : 0 will therefore evolve (approximately) into 1 : 1 : 1 at the detector, see, for example, [57, 58]. In Section 4.4, we will see that (4.1) can significantly change in the presence of new physics effects, which opens new possibilities to test such effects. However, (4.1) also implies that there could be some sensitivity to standard flavor mixing, which may be complementary to Earth-based experiments, see discussions in [59–71]. In the light of the current bounds for astrophysical neutrino fluxes from IceCube, however, such applications might be unlikely.

The main observable in neutrino telescopes is muon tracks from charged current interactions of muon neutrinos, producing Cherenkov light, which can be detected in so-called digital optical modules (DOMs). Because of the long muon range which is increasing with energy, the muon track does not have to be fully contained in the detector volume, which leads to excellent statistics increasing with energy. In addition, muon tracks have a very good directional resolution (order one degree). Additional event topologies include electromagnetic (mostly from electron neutrinos) and hadronic (from tau neutrinos) cascades, as well as neutral current cascades for all flavors. For even higher energies, the tau track may be separated, leading to so-called double bang or lollipop events; see [72, 73] for an overview. In practice, the main “flavor” analysis so far performed by the IceCube collaboration has been a cascade analysis [74]. To see that, consider that electromagnetic (from 𝜈𝑒) and hadronic (from 𝜈𝜏) cascades cannot be distinguished. A useful observable is therefore the ratio of muon tracks to cascades [59]̂‌𝑅≡𝜙𝜇𝜙𝑒+𝜙𝜏.(4.2)

Note that neutral current events will also produce cascades, and 𝜈𝜏 will also produce muon tracks in 17% of all cases, which, in practice, have to be included as backgrounds. In [74], the contribution of the different flavors to the cascade rate for a 𝐸−2𝜈 extragalactic test flux with equal contributions of all flavors at the Earth was given as electron neutrinos 40%, tau neutrinos 45%, and muon neutrinos 15% (after all cuts). This implies that charged current showers dominate and that electron and tau neutrinos are detected with comparable efficiencies, that is, that (4.2) is a good first approximation to discuss flavor at a neutrino telescope. The benefit of this flavor ratio is that the normalization of the source drops out. In addition, it represents the experimental flavor measurement with the simplest possible assumptions. For a pion beam source, one finds ̂‌𝑅≃0.5 at the detector, for a muon damped source ̂‌𝑅≃0.6, and for a neutron beam source ̂‌𝑅≃0.3, with some dependence on the mixing angles. In principle, this and other observables, such as the ratio between electromagnetic and hadronic showers, can therefore be used to determine the flavor composition of the sources, see, for example, [75–79] (unless there are new physics effects present, see for example, [80], which may yield similar ratios at the detector for different flavor compositions at the source).

We show the ratio ̂‌𝑅 as a function of energy for a GRB neutrino flux in Figure 7. At the “plateau” of the flux, where the largest number of events is expected, one can clearly see the flavor transition between a pion beam and a muon damped source at higher energies. In addition, the neutron decays have an impact at very low energies. In this figure, the uncertainty on ̂‌𝑅(3𝜎) coming from the current mixing angle uncertainties is shown, as well as the expectation for the next generation of reactor and long-baseline experiments (2015), and for a high-precision neutrino oscillation facility (2025). Obviously, the current uncertainties on the mixing angles are still too large to allow for a clear identification of the flavor composition at the source, whereas already the knowledge from the next generation will allow for a flavor ratio discrimination—at least in principle.

Figure 7 

Flavor ratio ̂‌𝑅 as a function of 𝐸𝜈 at the detector for a GRB neutrino flux with the Waxman-Bahcall (WB) “plateau” (in 𝐸2𝜙𝜇) between 105 and 107 GeV (cf. [17]). The shaded regions show the impact of the 3𝜎 mixing angle uncertainties now (current) [81], in about 2015 (next generation, dominated by Daya Bay and T2K [82]), and in about 2025 (for a neutrino factory [83, 84]). Here 𝜃13=0 is assumed for the sake of simplicity, figure taken from [30].

4.2. Detector Response and Impact of Spectral Shape

For time-integrated point source searches in IceCube [13], the simplest possible approach to describe the event rate of muon tracks is to use the exposure Exp(𝐸𝜈,𝛿)≡𝐴eff𝜈(𝐸𝜈,𝛿)𝑡exp, where 𝐴eff𝜈 is the neutrino effective area and 𝑡exp is the observation (exposure) time. Here 𝐴eff𝜈(𝐸𝜈,𝛿) is a function of the flavor or interaction type (which we do not show explicitly), the incident neutrino energy 𝐸, and the declination of the source 𝛿. The neutrino effective area already includes Earth attenuation effects (above PeV energies) and event selection cuts to reduce the backgrounds, which depend on the type of source considered, the declination, and the assumptions for the input neutrino flux, such as the spectral shape. Normally, the cuts are optimized for an 𝐸−2𝜈 flux. The total event rate of a neutrino telescope can be obtained by folding the input neutrino flux with the exposure as𝑁=∫𝑑𝐸𝜈Exp(𝐸𝜈,𝛿)𝜙𝜇(𝐸𝜈)=∫𝑑𝐸𝜈𝐴eff𝜈(𝐸𝜈,𝛿)𝑡exp𝜙𝜇(𝐸𝜈).(4.3)

Here 𝜙𝜇(𝐸𝜈) is, for point sources, given in units of GeV⁻¹ cm⁻² s⁻¹ for neutrinos and antineutrinos added. If backgrounds are negligible, the 90% (Feldman-Cousins) sensitivity limit 𝐾90 for an arbitrarily normalized input flux used in (4.3) can be estimated as 𝐾90∼2.44/𝑁 [85]. This implies that a predicted flux at the level of the sensitivity limit, irrespective of the spectral shape, would lead to the same number (2.44) of events. The 90% confidence level differential limit in terms of 𝐸2𝜈𝜙𝜇 can be defined as 2.3𝐸𝜈/Exp(𝐸𝜈,𝛿), see, for example, [86].

The interplay between spectral shape and detector response has, for instance, been studied in [46]. In order to illustrate that, we show in Figure 8 the limits for selected muon neutrino spectra and two different source declinations, corresponding to downgoing and upgoing events, for the time-integrated point source search in IC40. The limits are computed with (4.3), where the normalization is obtained from 𝐾90∼2.44/𝑁 and the spectral shape is predicted by the model. In this figure, the thick horizontal lines show the limits for an 𝐸−2𝜈 flux (in the dominant energy range), and the thick curves the differential limits. One can easily see that the differential limits are useful, since any “well behaved” neutrino flux will stay below these limits. In addition, the differential limit shows the energy range where the instrument is most sensitive to a flux, whereas, for the horizontal lines, representing an 𝐸−2𝜈 flux limit, only the contribution close to the differential limit minimum contributes to (4.3). In Figure 8, the minimal and maximal energies for the horizontal lines are indeed arbitrarily chosen, since the limit hardly depends on these (for the chosen example, this holds as long as the energy range around the differential limit minimum within about 1.5 orders of magnitude in energy is included. Sometimes the range where 90% of the events come from is shown). From the different fluxes in Figure 8 it is clear that the interplay between spectral shape and detector response is important. For example, spectrum no. 4 will be better constrained in Figure 8(a) (downgoing events) than in Figure 8(b) (upgoing events) because the differential limit peaks at higher energies—in spite of the lower absolute performance at the differential limit minimum. The reason is the coincidence between spectral peak and differential limit minimum. This picture changes completely in the presence of strong magnetic field effects on the secondaries, see spectrum no. 2, where even the imprint of these effects becomes important. The comparison to the 𝐸−2𝜈 test flux (horizontal line) clearly demonstrates that the 𝐸−2𝜈 assumption is insufficient for all practical purposes, since it does not take into account the energy range of the flux. In addition, as we discussed around (3.5) earlier, the 𝐸−2𝜈 assumption only holds for a very special case for the target photon density. Therefore, detector response and source model are intimately connected, and it is not quite clear if some sources may be missed just because the detector response does not match the source prediction. We will come back to this issue in Section 5.

(a)

(b)

Figure 8 

Limits for selected muon neutrino spectra (including flavor mixing) from [34] for different declination bands for the time-integrated point source search in IC40 (90% CL) [13]. The thick lines show the limits for an 𝐸−2𝜈 flux (in the dominant energy range), and the thick curves the differential limits, figure taken from [46]. The different spectra correspond to the test points in Figure 10 (same numbers used).

Figure 8 is also useful to illustrate the impact of the source declination, which is, for IceCube, a measure for the direction of the muon track in terms of the nadir angle because of the location at the South Pole. Note, however, that for other neutrino telescopes, such as ANTARES, this relationship is nontrivial. Obviously, the differential limit is different in the left and right panels of this figure, corresponding to downgoing and upgoing events, respectively. First of all, note that the downgoing events have to fight the atmospheric muon background, which leads to a worse performance at low energies because of appropriate cuts. However, the upgoing events suffer from the Earth attenuation for energies ≫PeV, which leads to the better sensitivity of the downgoing events for high energies, and the shift of the differential limit minimum to higher energies. From the discussion above it is clear that both event types (upgoing and downgoing) are complementary not only because they test a different part of the sky, but also because they test different energy ranges.

Finally, let us briefly comment on the statistics expected for neutrinos. First of all, it is clear from the discussion above that per definition any current IC40 limit, obtained over about one year of data taking, is compatible with about 2.4 events at the current 90% confidence limit. Assuming that the effective area increases by about a factor of three to four from IC40 to IC86 [87], one can extrapolate that the current limits are compatible with 2.4⋅4⋅10≃100 events over ten years of full IceCube operation if the current bounds are saturated. Any further nonobservation of events will reduce this maximal expectation. Therefore, one can already say right now that any conclusion about the astrophysical neutrino sources, the sources of the cosmic rays, or leptonic versus hadronic models for 𝛾-ray observations will most likely be based on the information from many sources of one class. A typical example is the stacking of GRBs using their gamma-ray counterparts, such as in [15]. The aggregation of fluxes, no matter if diffused or stacked, will however imply new systematics and model-dependent ingredients, see discussion in [30].

4.3. Glashow Resonance to Discriminate 𝑝𝑝 from 𝑝𝛾?

A useful observable may be the Glashow resonance 𝜈𝑒+𝑒−→𝑊−→anything at around 6.3 PeV [34, 57, 88–92] to distinguish between neutrinos and antineutrinos in the detector, since this process is only sensitive to 𝜈𝑒. For photohadronic (𝑝𝛾) interactions, however, mostly 𝜋+ and therefore 𝜈𝑒 are produced at the source, see (2.1) and (2.2), which means that no excess of events should be seen at this resonance energy—at least in the absence of flavor mixing. On the other hand, for 𝑝𝑝 interactions in the source, 𝜋+ and 𝜋− are produced in about equal ratios, which increases the 𝜈𝑒 production rate. Therefore, the Glashow resonance is frequently proposed as a discriminator between 𝑝𝑝 and 𝑝𝛾 sources. Note that this argument is especially interesting for 𝐸−2𝜈 neutrino fluxes, whereas any other, (significantly different) spectral index may be a clear indicator for 𝑝𝛾 interactions.

In [34], the electron neutrino-antineutrino ratio at source (Section 3.3) and detector (Section 4.3) has been computed explicitly. At the source, the following observations have been made.(i) Additional production modes (such as direct and multipion production) in photohadronic interactions produce 𝜋− in addition to 𝜋+. The pure 𝜋+ source therefore does not exist, and up to 20% contamination from 𝜈𝑒 at the source has to be accepted even in this case. As a consequence, the Glashow resonance must be seen to some degree, even for the 𝑝𝛾 source.(ii) Since the photohadronic interaction in (2.1) also exists for neutrons producing 𝜋− in this case, any optical thickness to neutron escape will lead to a 𝜋+-𝜋− equilibration. This means that the optically thick 𝑝𝛾 source cannot be distinguished from the 𝑝𝑝 source by the Glashow resonance.(iii) Neutron decays, which are inherently present in any photohadronic source, lead to 𝜈𝑒 flux, faking a 𝜋− contribution and therefore a 𝑝𝑝 source in particular energy intervals (determined by the maximal proton energy).

As a consequence, only the “𝑝𝛾 optically thin” (to neutron escape) source might be uniquely identified if less than 20% of 𝜈𝑒 contamination are found. On the other hand, one cannot uniquely identify a 𝑝𝑝 source. Note that modern approaches take the composition between 𝑝𝑝 and 𝑝𝛾 interactions as a variable, see [90–92]. However, these approaches cannot typically describe the contamination from neutron decays, because this depends very much on the model.

At the detector, the use of the Glashow resonance becomes even more complicated for the following reasons.(i) Flavor mixing reintroduces a 𝜈𝑒 component from 𝜈𝜇 produced in 𝜇+ decays even for pure 𝜋+ production, compare (2.2).(ii) As a consequence, the flavor composition at the source is important and has to be determined at the same time. For example, a muon damped 𝑝𝑝 source may be easily mixed up with a pion beam 𝑝𝛾 source.(iii) The Glashow resonance occurs at a specific energy (6.3 PeV), which means that for this process the energy dependence is not important and only one particular energy matters. However, the transformation of the energy from source to detector depends on redshift and a possible Lorentz boost, compare (3.11), which have to be known to draw any conclusions.

In summary, the discovery of a neutrino signal with significantly suppressed Glashow resonance events will be interesting and may allow for possible conclusions about the source if flavor composition, 𝑧, and Γ are known. On the other hand, the detection of a pronounced Glashow resonance events is probably nonconclusive for the physics of the source.

4.4. Testing New Physics in the Neutrino Propagation

In the Standard Model, the transition probability 𝑃𝛼𝛽 in (3.11) is described by the usual flavor mixing in (4.1), which is independent of energy. New physics effects may lead to deviations from this picture, where the effects discussed in the literature include sterile neutrinos, neutrino decay, quantum decoherence, Lorentz invariance violation, among others, see [53] for a review. In some of these cases, 𝑃𝛼𝛽 may even be a function of energy 𝑃𝛼𝛽(𝐸). We choose neutrino decay as an example in this section, see [38, 80, 90, 93–100], where in [98–100] energy-dependent effects have been considered.

Following [100], for decay into invisible particles (see, e.g., [90, 101] for the systematical discussion/classification of possible scenarios. In general, neutrinos may decay into particles invisible for the detector, or other visible active flavors. Furthermore, the decay may be complete, i.e., all particles have decayed, or incomplete, where the energy dependence of the decay can be seen. For the case of complete decays, there are only eight different effective decay scenarios (each mass eigenstate can be stable or unstable), no matter if the decays are visible or invisible or if intermediate (unstable) invisible states are involved [90]. Likewise, there are only eight scenarios for invisible incomplete decays, whereas the treatment of incomplete visible decays is more complicated, see, e.g., [101, 102]), the transition probability can be described by a modified version of (4.1):𝑃𝛼𝛽=3∑𝑖=1|𝑈𝛽𝑖|2|𝑈𝛼𝑖|2𝐷𝑖(𝐸)with𝐷𝑖(𝐸)=exp(−̂𝛼𝑖𝐿𝐸),(4.4) as the damping coefficient [103]. Here ̂𝛼𝑖=𝑚𝑖/𝜏0𝑖 with 𝜏0𝑖 the rest frame lifetime for mass eigenstate 𝜈𝑖. Typically the neutrino lifetime is quoted as 𝜏0𝑖/𝑚𝑖 since 𝑚𝑖 is unknown, see [53] for a review. Note that this transition probability is energy dependent, compared to (4.1). In this case, the flavor ratio ̂‌𝑅 in (4.2) can be rewritten aŝ‌𝑅=𝑃𝑒𝜇(𝐸)̂‌𝑋(𝐸)+𝑃𝜇𝜇(𝐸)[𝑃𝑒𝑒(𝐸)+𝑃𝑒𝜏(𝐸)]̂‌𝑋(𝐸)+[𝑃𝜇𝑒(𝐸)+𝑃𝜇𝜏(𝐸)](4.5) if ̂‌𝑋(𝐸)≡𝑄𝑒(𝐸)/𝑄𝜇(𝐸) is the ratio between electron neutrinos and muon neutrinos ejected at the source (assuming that hardly any tau neutrinos are produced). Equation (4.5) carries now two energy dependencies: the energy dependence of the flavor composition at the source ̂‌𝑋(𝐸), and the energy dependence of the new physics effect 𝑃𝛼𝛽(𝐸), which have to be intrinsically disentangled. On the other hand, the energy dependence of the new physics effect may provide a unique signature [98, 99], and the energy dependence of the flavor composition at the source may help to disentangle different scenarios [100].

In general, there are 23=8 decay scenarios for invisible incomplete decays, since every (active) mass eigenstate may be either stable or unstable. In Figure 9(a), ̂‌𝑅 is shown as a function of the initial flavor composition ̂‌𝑋 for these different (complete) decay scenarios, where filled disks correspond to stable mass eigenstates and unfilled disks to unstable mass eigenstates (the scenario with only unstable states is of course not shown, since no neutrinos can be detected in this case). In this panel, the different types of sources are also marked. One can clearly see that especially for a pion beam source, three of the scenarios cannot be distinguished, whereas these can be disentangled in principle for any other type of source. In addition, for three scenarios (the scenarios where only one mass eigenstate is stable), ̂‌𝑅 is independent of the initial flavor composition. Note that the scenario with the largest ̂‌𝑅 (𝑚1 and 𝑚2 unstable) faces the strongest constraints because of the observation of neutrinos from supernova 1987A.

(a)

(b)

Figure 9 

(a) ̂‌𝑅 as a function of the initial flavor composition ̂‌𝑋 for different decay scenarios (complete decays), where filled disks correspond to stable mass eigenstates and unfilled disks to unstable mass eigenstates. (b) ̂‌𝑅 as a function of energy for a pion beam to muon damped source (cf. Figure 5, TP 13) for a specific decay parameter (same for all unstable states), as given in the panel. The dashed curves stand for complete decays, figures adopted from [100].

In Figure 9(b), ̂‌𝑅 is shown as a function of energy for a particular choice of the decay parameter times distance and a specific source with a flavor transition from pion beam to muon damped source at about 106 GeV. First of all note that because of the exponential damping in (4.4), the decays are practically complete for energies 𝐸≪108 GeV, and the neutrinos are stable for energies 𝐸≫108 GeV for the chosen decay parameters. This is why the scenarios deviate from the complete decay curves (dashed) at very high energies and converge into the standard scenario there (all neutrinos stable). The pion beam (low energies) cannot distinguish three of the scenarios at low energies, as expected. However, above 106 GeV, where the flavor transition occurs into a muon damped source, the corresponding curves split up (and the neutrino decays are still complete), before they converge into the standard case. This example illustrates that the energy-dependent flavor transition of a specific source may be a useful for new physics tests, provided that enough statistics can be collected.

5. Application to Generic (AGN-Like) Sources

As it is illustrated in Figure 6 (see also (3.1)), the proton 𝑁′𝑝 and photon densities 𝑁′𝛾 within the source control the secondary meson, and hence the neutrino production. In this section, we follow the ansatz in [34]: we assume that the target photons are produced in a self-consistent way, by the synchrotron emission of coaccelerated electrons (positrons); see, for example, [104] for a corresponding specific (BL Lac) AGN blazar model. The main purpose of this model is the prediction of spectral shape and flavor composition of the neutrino source, while it cannot predict the normalization of the neutrino flux in the form presented here. In addition, no neutrinos have been observed yet, which means that an ansatz tailor-made for neutrinos may be useful to test the interplay between detector response and source model; after all, there may be sources for which the optical counterpart is absorbed, so-called “hidden sources,” see, for example, [105–108]. The ingredients to this model are comparable to the conventional GRB neutrino models, discussed in Section 6, where however the origin of the target photons is different.

5.1. Additional Model Ingredients

The primaries, in this case protons and electrons, are assumed to be injected with an (𝐸′)−𝛼 injection spectrum, where the (universal) injection index 𝛼 is one of the important model parameters. They are assumed to lose energy dominantly by synchrotron radiation, controlled by 𝐵′, and adiabatic cooling, controlled by 𝑅′, which also determine their maximal energies. This means that (3.9) is applied to the primaries, which implies that the acceleration and radiation zones are different in this model. As a consequence, the electron spectrum becomes loss-steepened by one power.

The synchrotron photons (the target photons in (3.1)) are computed in the Melrose-approximation [109] averaged over the pitch angle. The power radiated per photon energy 𝜀′ by one particle with energy 𝐸′, mass 𝑚 and charge 𝑞 in a magnetic field 𝐵′ is given by𝑃′𝑠(𝜀′,𝐸′)=1.8⋅√3𝑞3𝐵′16𝜀0𝑚𝑐ℎ⋅(𝜀′𝜀′𝑐)1/3𝑒−𝜀′/𝜀′𝑐with𝜀′𝑐=3𝑞𝐵′16𝑚(𝐸′𝑚𝑐2)2.(5.1)

We have to convolute this with the spectrum of radiating electrons as𝑃′(𝜀′)=∫∞0𝑑𝐸′𝑁′𝑒(𝐸′)𝑃′𝑠(𝜀′,𝐸′).(5.2)

The number of produced photons per time can be computed with𝑄′𝛾(𝜀′)=𝑃′(𝜀′)𝜀′.(5.3)

Assuming that the photons escape (and hardly interact again), the steady photon spectrum, which is needed for the computation of photohadronic interactions, can be estimated by multiplying 𝑄′𝛾 with escape time 𝑡′esc≃𝑅′/𝑐 of the photons. The synchrotron spectral index (𝜀′)−𝛼𝛾 obtained from this approach is 𝛼𝛾=(𝛼𝑒−1)/2+1=𝛼/2+1, which means that the dependence on the primary injection index 𝛼 is small. The rest of the computation follows Section 3.

The main parameters in this model are 𝛼, 𝑅′, and 𝐵′. For the sake of simplicity, we assume in the following that the source is only moderately Lorentz boosted with respect to the observer’s frame, that is, 𝑅′≃𝑅 and 𝐵′≃𝐵. In this case, a convenient description of the parameter space of interest is the Hillas plot [110]. In order to confine a particle in a magnetic field at the source, the Larmor radius has to be smaller than the extension of the acceleration region 𝑅. This can be translated into the Hillas condition for the maximal particle energy𝐸max[GeV]≃0.03⋅𝜂⋅𝑍⋅𝑅[km]⋅𝐵[G].(5.4)

Here 𝑍 is the charge (number of unit charges) of the accelerated particle, 𝐵 is the magnetic field in Gauss, and 𝜂 can be interpreted as an efficiency factor or linked to the characteristic velocity of the scattering centers. Potential cosmic ray sources are then often shown in a plot as a function of 𝑅 and 𝐵, as it is illustrated in Figure 10(a), by the numbered disks (see legend for possible source correspondences). Assuming that a source produces the highest energetic cosmic rays with 𝐸≃1020 eV, one can interpret (5.4) as a necessary condition excluding the region below the Γ=1 line in Figure 10 (for protons with 𝜂=1). If one allows for relativistic boosts between source and observer, this condition is relaxed, as one can read off from the figure (in this case, 𝑅 and 𝐵 in the plot have to interpreted as 𝑅′ and 𝐵′). However, this method does not take into account energy loss mechanisms, which lead to a qualitatively different picture, see, for example, [111, 112], and which are implied in our model. In the following, we will study the complete parameter space covered by Figure 10(a) without any prejudice. Since the location of the sources in Figure 10 cannot be taken for granted, we will refer to the individual sources as “test points” (TP), and leave the actual interpretation to the reader.

(a)

(b)

Figure 10 

(a) Possible acceleration sites in Hillas plot as a function of 𝑅 and 𝐵 (version adopted from M. Boratav). (b) Classification of sources for injection index 𝛼=2 in this plot (see main text). Some points from left plot are shown for orientation, as well as two new points (12 and 13), figure taken from [34].

Concerning the limitations of the model, it certainly does not apply exactly to all types of sources. For example, in supernova remnants, 𝑝𝑝 (proton-proton) or 𝑝𝐴 (proton-nucleus) interactions may dominate the neutrino production, which would require additional parameters to describe the target protons or nucleons. In addition, at ultra-high energies, heavier nuclei may be accelerated. The spirit of this model is different: it is developed as the simplest (minimal) possibility including nontrivial magnetic field and flavor effects. Another ingredient is the target photon density, which is assumed to come from synchrotron emission of coaccelerated electrons here. In more realistic models, typically a combination of different radiation processes is at work. However, in many examples with strong magnetic fields, the specific shape of the photon spectrum is less important for the neutrino spectral shape than the cooling and decay of the secondaries, which depend on particle physics only. To check this, we have tested the hypothesis that acceleration and radiation zones of the electrons are identical (i.e., the electron spectrum is not loss-steepened by synchrotron losses, which is actually a simpler version of this model). Thus, while it is unlikely that the model applies exactly to a particular source, it may be used as a good starting hypothesis.

5.2. Flavor Composition at the Source

We discussed the flavor composition at the source already in Section 3.3, where we also showed several examples for the energy dependence in Figure 5. Let us now approach this from a systematical point of view as a function of 𝑅 and 𝐵 (for 𝛼=2). In this case, a qualitative classification of the sources can be found in Figure 10(b), where it is implied that a certain flavor ratio can be clearly identified over one order of magnitude in energy close enough to the spectral peak; see [34]. From this figure, we find that for 𝐵≳1012 G, all charged species lose energy so rapidly that the neutron decays dominate the neutrino flux. For 𝐵≳1010 G, several processes compete, leading to an undefined source. For 𝐵≲1 kG, the sources behave as classical pion beams, which typically applies to sources on galactic scales. In the intermediate range, 1kG≲𝐵≲1010G, the source classification somewhat depends on the spectral shape, since this affects possible muon pile-up effects, the energy range close to the spectral peak, and the competition of several effects. Depending on 𝑅, muon beam to muon damped sources (as a function of energy), muon damped sources, and pion beam to muon damped sources are found. There is some dependence on 𝛼, which affects the spectral shapes. For instance, for 𝛼∼4 a similar pion spectrum as for GRBs is obtained, which leads to a pion beam to muon damped sources for typical GRB parameters in consistency with [37, 38]. In summary, the pion beam source assumption is safe if 𝐵≲1 kG, whereas for stronger magnetic fields the magnetic field effects on the secondaries have to be taken into account. However, for given parameters, this flavor composition can be predicted. For instance, 𝑅 may be estimated from the time variability of the source, 𝐵 from energy equipartition arguments, and 𝛼 from the observed photon spectrum.

5.3. Interplay between Spectral Shape and Detector Response

We discussed already earlier in Section 4.2 the interplay between spectral shape and detector response, see [46]. From Figure 8, it is clear that spectral shapes with a peak at the position of the differential limit minimum can be better limited than others. This is especially clear if two event types or detectors are compared for the same spectrum. However, how should one quantify that for different spectra seen by the same detector? Consider, for example, the fluxes no. 2 and no. 11 in Figure 8(a), both leading to the same event rate by definition. Which of the two neutrino sources is the detector more sensitive to in terms of the physics of the source? In order to address this question, it is useful to assign a single number to each spectrum which measures how much energy in neutrinos can be tested for a specific spectrum and event type. We choose the energy flux density𝐹𝛽=∫𝐸𝜈𝜙𝛽𝑑𝐸𝜈,(5.5) as this quantity, which we show in units of erg cm⁻² s⁻¹ for point sources in order to distinguish it from 𝐸2𝜈𝜙 in units of GeV cm⁻² s⁻¹ (1erg≃624 GeV).

This quantity measures the total energy flux in neutrinos, and it is useful as a performance indicator measuring the efficiency of neutrino production in the source. In order to see that, consider the alternative derivation of 𝐹𝛽 from (3.11) (neglecting a possible Lorentz boost of the source)𝐹𝛽=∑𝛼=𝑒,𝜇,𝜏𝐿𝜈𝛼4𝜋𝑑2𝐿,where𝐿𝜈𝛼=𝑉′∫𝐸′𝜈𝑄′𝜈𝛼𝑑𝐸′𝜈,(5.6) where is the “neutrino luminosity” and 𝑉′ is the volume of the interaction region. Since the neutrinos originate mostly from pion decays and take a certain fraction of the pion energy (about 1/4 per produced neutrino for each charged pion), the neutrino luminosity is directly proportional to the (internal) luminosity of protons 𝐿int (or the proton energy dissipated within a certain time frame Δ𝑇) and the fraction of the proton energy going into pion production, commonly denoted by 𝑓𝜋 (if the energy losses of the secondaries can be neglected). Since a possibly emitted photon flux can be often linked to 𝐿int by energy partition arguments, one has 𝐹𝛽∝𝑓𝜋×𝐿int∝𝑓𝜋×𝐿𝛾, and 𝐹𝛽 is a measure for 𝑓𝜋×𝐿int of the source (if no photon counterpart is observed), or even 𝑓𝜋 itself (if a photon counterpart is observed).

Regions for a specific sensitivity to 𝐹𝛽 are shown as a function of 𝑅 and 𝐵 in Figure 11 for 𝜈𝜇(Figure 11(a)) and 𝜈𝜏(Figure 11(b)), for several source declinations in IceCube and for Earth-skimming tau neutrinos in Auger (in this case, for a diffuse flux). There are several conclusions from this figure: first of all, IceCube responds very well to the usual suspects, such as AGNs (Figure 11(a), points 3 to 6). Even for 𝜈𝜏, which produce a muon track in only 17% of all cases, most of this parameter space can be covered (Figure 11(b)). Thus, it is clear that most sources will be also detected if the partition between 𝜈𝜇 and 𝜈𝜏 was heavily disturbed, such as by a new physics effect in the neutrino propagation. The sensitive region, however, somewhat depends on the source declination. For very large values of 𝐵, the neutrino energies are lower and instruments such as the DeepCore array may respond better to the flux (Figure 11). On the other hand, the region where the UHECRs are expected in this model (lower right corner) is better covered by Auger (Figure 11(b)). This is not a big surprise: the neutrino spectrum follows the proton spectrum in the absence of strong magnetic field effects on the secondaries, which means that the spectra (cf. Figure 8) peak at high energies. A very interesting region may be the gap between the IceCube and Auger regions: perhaps a future instrument such as KM3NeT should optimize their geometry to be complementary in terms of energy range coverage.

(a)

(b)

Figure 11 

Regions where the sensitivity in 𝐹𝛽 exceeds 10−9 erg cm−2 s−1 [sr−1] for 𝛽=𝜇 (a) and 𝛽=𝜏 (b) for several selected data samples (90% CL, 𝛼=2). The dashed regions “UHECR” indicate where 1020 eV cosmic ray protons are expected to be produced in the model. The different sources (test points) correspond to the ones in Figure 10, figure taken from [46].

6. GRB Neutrinos and the Multimessenger Connection

Recall again that, as it is illustrated in Figure 6 (cf. (3.1)), the proton 𝑁′𝑝 and photon densities 𝑁′𝛾 control neutrino production. In the previous section, we assumed that 𝑁′𝛾 is produced by synchrotron emission of coaccelerated electrons. Here we emphasize the multimessenger connection, assuming that 𝑁′𝛾 can be derived from the gamma-ray observation [30, 44]; see also [15, 20, 21, 113] for analytical approaches. The motivation is also different from the previous section: whereas we were interested in a systematic parameter space study of spectral shape and flavor ratio before, the main emphasis here is the prediction of spectral shape, flavor composition, and absolute neutrino flux normalization for a specific set of GRBs observed in gamma rays.

6.1. Additional Model Ingredients

Here we describe the key ingredients of the conventional (numerical) fireball GRB model for neutrino emission, following [30], where we focus on the normalization. For a detailed comparison to the analytical calculations, see [29, 44]. Because the model does not describe the neutrino production in a time-resolved way, it makes sense to relate the neutrino production to the (bolometric) gamma-ray fluence 𝑆bol (in units of erg cm⁻²) during the burst. The isotropic bolometric equivalent energy 𝐸iso,bol (in erg) in the source (engine) frame can be then obtained as𝐸iso,bol=4𝜋𝑑2𝐿(1+𝑧)𝑆bol.(6.1)

It can easily be boosted into the SRF by 𝐸′iso,bol=𝐸iso,bol/Γ. Assuming energy equipartition between photons and electrons, the photons carry a fraction 𝜖𝑒 (fraction of energy in electrons assuming equipartition between kinetic energy in electrons and energy in photons) of the total energy 𝐸iso,tot as𝐸iso,tot=𝜖−1𝑒⋅𝐸iso,bol.(6.2)

In order to compute the photon and proton densities in the SRF, it turns out to be useful to define an “isotropic volume” 𝑉′iso𝑉′iso≃4𝜋(2Γ2𝑐𝑡𝑣(1+𝑧)⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟𝑅𝐶)2⋅(Γ⋅𝑐(𝑡𝑣(1+𝑧))⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟Δ𝑑′)∝Γ5,(6.3) where 𝑅𝐶 is the collision radius, Δ𝑑′ is the shell thickness of the colliding shells, and 𝑡𝑣 is the (observed) variability timescale. Thus, 𝑉′iso can be estimated from 𝑡𝑣, 𝑧, and Γ and can be regarded as the volume of the interaction region assuming isotropic emission by the source. Because of the intermittent nature of GRBs, the total fluence is assumed to be coming from 𝑁≃𝑇90/𝑡𝑣 such interaction regions, where 𝑇90 is the duration of the burst (time during which 90% of the total energy is observed). Now one can determine the normalization of the photon spectrum in (3.1) from∫𝜀′𝑁′𝛾(𝜀′)𝑑𝜀′=𝐸′iso,bol𝑁𝑉′iso,(6.4) assuming that the spectral shape is determined by the observed spectrum. Similarly, one can compute the normalization of the proton spectrum in (3.1) by∫𝐸′𝑝𝑁′𝑝(𝐸′𝑝)𝑑𝐸′𝑝=1𝑓𝑒𝐸′iso,bol𝑁𝑉′iso,(6.5) where 𝑓𝑒 is the ratio between energy in electrons and protons (𝑓−1𝑒: baryonic loading). Note that, in the end, we will obtain the neutrino flux 𝜙 per time frame per interaction region from (3.11) with ̂‌𝑁=𝑉′iso. Assuming that the magnetic field carries a fraction 𝜖𝐵 of 𝐸iso,tot, one has in addition𝑈′𝐵=𝜖𝐵𝜖𝑒⋅𝐸′iso,bol𝑁𝑉′isoor𝐵′=√8𝜋𝜖𝐵𝜖𝑒⋅𝐸′iso,bol𝑁𝑉′iso.(6.6)

Typical values used in the literature are 𝑓𝑒∼𝜖𝑒∼𝜖𝐵≃0.1 (see, e.g., [21]). An explicit calculation for 𝐵′ yields𝐵′≃220(𝜖𝐵𝜖𝑒)1/2(𝐸iso,bol1053erg)1/2(Γ102.5)−3×(𝑡𝑣0.01s)−1(𝑇9010s)−1/2(1+𝑧3)3/2kG.(6.7)

In summary, once the photon (from observation) and proton (typically 𝐸−2𝑝) spectral shapes are determined, the proton and photon densities in the source and 𝐵′ can be calculated with the above formulas from the observables (gamma-ray fluence, Γ, 𝑧, 𝑡𝑣). Equation (6.3) implies that for fixed 𝑡𝑣, the larger Γ, the larger the interaction region, and the smaller the photon density in (6.4), which directly enters the fraction of fireball proton energy lost into pion production 𝑓𝜋∝Γ−4 [20], or, consequently, 𝐸2𝜈𝜙∝Γ−4. Therefore, the main contribution to the neutrino flux is often believed to come from bursts with small Lorentz factors, see discussions in [20, 114, 115]. We follow this conventional fireball approach in the following. Note that this numerical approach contains all the ingredients of [20, 21] explicitly—such as the cooling of the secondaries. The neutrino emission can be then easily computed as shown in Figure 6. In addition, note that there are alternatives to the model. For instance, if the bursts are alike in the comoving frame, as suggested in [116], one has 𝐸2𝜈𝜙∝Γ2 [30]. One can read off from (6.3) that the correlation Γ∝𝑡−3/5𝑣 is expected in that case, since 𝑉′iso will be similar for the bursts. See [30] for a discussion of the neutrino flux for different model hypotheses, and how the neutrino flux can in principle be used to discriminate among these.

6.2. Systematics in the Interpretation of Aggregated Fluxes

Since the number of neutrinos expected from a single GRB is small, dedicated aggregation methods are needed. For instance, one may search for the diffuse flux from GRBs, which, however, has to fight the background from atmospheric neutrinos. Another possibility is to use the gamma-ray observation to infer on time window and direction of the neutrino signal, which effectively leads to significantly reduced backgrounds. In addition, a procedure such as the one in Section 6.1 may be used to predict the absolute neutrino flux, its shape (see, e.g., Figure 4), or its flavor composition (see, e.g., Figure 7) on a burst-by-burst basis. Summing over many observed bursts, such an analysis is also called stacking analysis, see, for example, [15] for a recent example. The diffuse limit can be extrapolated from such a stacked flux, which is also called the quasi-diffuse limit; see [7] for details. Here we discuss some of the implications when a stacked neutrino flux is translated into a quasi-diffuse flux.

In order to illustrate this problem, following [30], let us consider the redshift distribution of the GRBs as an example. We show in Figure 12(a), a population of 10 000 bursts, representative for the number of GRBs in the visible universe over about 10 years (cf. histogram). These bursts are assumed to follow the star formation rate from Hopkins and Beacom [117] with the correction ℰ(𝑧) from Kistler et al. [118]. For the sake of simplicity, assume that all bursts have the same isotropic luminosity. From the discussion after (3.11), we have 𝐸2𝜈𝜙∝1/𝑑𝐿(𝑧)2, which means that closer GRBs, which are however rarer, will lead to a larger neutrino flux. Thus, it is the product 𝑑−2𝐿𝑑̇‌𝑁/𝑑𝑧 which determines the main contribution to the neutrino flux, shown as solid curve in Figure 12. While the peak contribution in terms of the GRB distribution is at z ~ 2–4, this contribution function peaks at 𝑧≃1. This observation has several implications: first of all, if the redshift 𝑧 is not measured, the neutrino flux may be overestimated if 𝑧=2 is assumed; compare (6.1). Second, the number of bursts contributing in the region 𝑧≃1 is rather small, which means that large statistical fluctuations are expected in quasi-diffuse flux estimates based on the stacking of a few bursts only.

(a)

(b)

Figure 12 

(a) Distribution of 10 000 bursts 𝑑̇‌𝑁/𝑑𝑧 as a function of redshift (histograms) and relative contribution of the individual GRBs 𝑑−2𝐿𝑑̇‌𝑁/𝑑𝑧 (solid curves). The dashed curve shows the exact distribution function. Here it is assumed that the GRBs follow the star formation rate from Hopkins and Beacom [117] with the correction ℰ(𝑧) from Kistler et al. [118]. (b) Probability that the quasi-diffuse flux extrapolated from a low statistics sample with 𝑛 bursts is larger than a certain fraction of the diffuse limit (see legend for different values of 𝑛). This function corresponds to (one minus) the cumulative distribution function of the probability density, figure taken from [30].

We quantify this systematical errors in Figure 12(b). In this panel, we show the probability that the quasi-diffuse flux extrapolated from a low statistics sample with 𝑛 bursts is larger than a certain fraction of the diffuse limit (see legend for different values of 𝑛). This function corresponds to (one minus) the cumulative distribution function of the probability density, and the step function corresponds to the diffuse limit. One can read off from this figure that for 𝑛≃100, corresponding to the analysis in [15], the quasi-diffuse extrapolation will be within 50% of the diffuse limit in the probability range corresponding to 90% of all cases (between 0.05 and 0.95). This means that a 50% error on the quasi-diffuse flux can be estimated from the redshift distribution only, while additional parameter variations increase this error [30].

6.3. Neutrino Flux Predictions from Gamma-Ray Observations

Since IceCube has not observed any GRB neutrino flux yet, there has been increasing tension between the model predictions [9, 20, 21] and the observation [15, 19] (recently, the superluminal propagation of neutrinos has also been proposed as a reason why no neutrinos have been seen, see, e.g., [119]). On the other hand, a direct comparison with the photomeson production in [9] was performed in [29] (see also Figure 1), demonstrating that the neutrino flux is actually underestimated in the analytical approaches. Therefore, the differences between the numerical approach in Section 6.1 and the analytical models in [20, 21] (based on [9]) have been identified in [44] by a recomputation of the analytical models and the analytical computation of a simplified version of the numerical code. As far as the astrophysics ingredients are concerned, these approaches can be shown to be equivalent, based on the same logic; see Section 6.1. The main differences are that magnetic field and flavor-dependent effects are explicitly included in numerical approach, additional pion, neutron, and kaon production modes are computed, and the full energy dependencies are taken into account.

We illustrate in Figure 13 the comparison between analytical and numerical approaches, where three different recent Fermi-measured GRBs have been chosen as examples: GRB 080916C, GRB 090902B, and GRB 091024. GRB 080916C has been selected, because it is one of the brightest bursts ever seen, although at a large redshift, and one of the best studied Fermi-LAT bursts. The gamma-ray spectrum of GRB 090902B can be fit by a Band function and a cutoff power law (CPL), which means that it can be used to illustrate the difference. GRB 091024 can be regarded as a typical example representative for many Fermi-GBM bursts [120], except from the long duration. Note that the first two bursts have an exceptionally large Γ≳1000, whereas Γ≃200 for the third burst. All three bursts have in common that the required parameters for the neutrino flux computation can be taken from the literature, in particular, the properties of the gamma-ray spectrum (including fluence), Γ, 𝑡𝑣, 𝑧, and 𝑇90, see figure caption for the references.

Figure 13 

Expected (time-integrated) neutrino flux for three different GRBs (parameters: see [120, 123] for GRB 080916C, [120, 124] for GRB 090902B, and [120, 125] for GRB 091024). First row: revision of IceCube analytical method conventional fireball calculation “CFB” (correction of shape, normalization from pion production efficiency 𝑐𝑓𝜋, and normalization from neutrino versus proton spectral shape 𝑐𝑆), leading to the revised fireball calculation “RFB.” Second row: comparison analytical (CFB, RFB) methods with simplified numerical method “WB Δ-approximation” (3.4) and full 𝑝𝛾 interactions. Third row: impact of adiabatic cooling on protons and secondaries. Fourth and fifth row: corresponding inverse timescales (rates) for protons and secondaries, respectively, courtesy of Hümmer [126].

In Figure 13, we show in the first row the computation of the predicted neutrino flux with the IceCube analytical method [21], called CFB (conventional fireball calculation) here. As described in [44], the corrections of shape, normalization from pion production efficiency 𝑐𝑓𝜋, and normalization from neutrino versus proton spectral shape 𝑐𝑆 (see also [45]) lead to a revised analytical calculation RFB. From the figure, one can easily read off that these revisions strongly depend on the burst parameters, especially the photon spectral shape. Comparing the predicted CFB fluxes with the bursts used for the IC40 analysis (Figure 1 in [15]), one can easily see that the expected fluxes of the first two bursts are about a factor 50 below that of the most luminous bursts in that analysis, and the third example about a factor of 5 below. This is expected from the scaling of the pion production efficiency ∝Γ−4 in that approach.

In the second row of Figure 13, we show the comparison of the analytical (CFB, RFB) methods with a simplified numerical method “WB Δ-approximation”, compare (3.4), and the full 𝑝𝛾 interactions. In most cases, the simplified numerical approach matches the method RFB rather well, which proves the validity of the derived corrections. However, for GRB 090902B (middle panel), the spectrum below the first break is different from the analytical estimate because the scalings of the weak decays limit the steepness of the spectrum there. The final numerical calculation including all production modes is then significantly enhanced again, especially due to multipion production [29]. Using a cutoff power law for GRB 090902B for the gamma-ray spectral fit (middle panel), the normalization of the prediction slightly reduces in that example because the photon density above the photon break is suppressed. Comparing the original CFB method with the final numerical computation “full 𝑝𝛾”, it is interesting that this can significantly deviate in both normalization and shape, but this deviation depends on the burst parameters. Very interestingly, a similar neutrino flux normalization for all three bursts is obtained, which means that it is probably not warranted to say that the neutrino flux from high-Γ bursts such as GRB 080916C is expected to be small. However, note that such extreme bursts only make up for a small fraction of the observed bursts, and the conclusions from neutrinos will be determined by the statistical properties of the burst sample in the stacking analysis. Finally, note that the numerical calculations in Figure 13 do not depend on any approximations, whereas different analytical methods lead to different predictions, similar to CFB. Therefore, the numerical computations should be regarded as the benchmark which defines the corrections, not vice versa. Within the simplest fireball neutrino model, there are only small model dependencies within the numerical approach. For instance, the integral limits in (6.5) (the minimal and maximal proton energies) have to be specified (note, however, that for an (𝐸′𝑝)−2 injection spectrum, the energy partition only logarithmically depends on the minimal and maximal proton energies), and a bolometric correction may have to be applied to (6.4)—which typically has small effects.

One of the effects included in the final result of the second row of Figure 13 is the adiabatic cooling of primaries (protons) and secondaries (muons, pions, kaons). We illustrate the effect of this cooling component in the lower three rows of the figure, where we show the impact on the final result (third row) and the respective inverse timescales (rates) for protons (fourth row) and secondaries (fifth row). For the protons (fourth row), synchrotron losses are assumed to determine the maximal energies in the absence of adiabatic cooling, whereas the larger of the synchrotron or adiabatic cooling loss rates determines the maximal proton energy otherwise. From the fourth row, one can also easily read off that the energy losses due to 𝑝𝛾 interactions are typically subdominant. The comparison with the third row illustrates that, depending on the burst parameters, the kaon hump (the rightmost one) can be suppressed by the adiabatic cooling of the protons, whereas the normalization of the spectra is hardly affected. From the fifth row, one can also read off that adiabatic cooling may have a small effect on the muons, for which it sometimes dominates in a small energy range. This leads to a small suppression of the first hump, coming from the muon decays (see left and middle panels of third row).

We show in Figure 14 the predicted quasi-diffuse neutrino flux from the above numerical method applied to the IC40 bursts for the same bursts and parameters (solid black curve) used in that analysis, which is about one order of magnitude lower than the IC40 limit and a factor of two below the current limit. In this figure a number of systematical errors are shown as well, such as the statistical error discussed in the previous subsection, and the estimated astrophysical uncertainty (by varying the unknown parameters, such as proton injection index 𝛼𝑝=1.8⋯2.2, variability timescale 𝑡𝑣 by one order of magnitude around the IceCube standard values, 𝑡𝑣=0.01 s for long bursts, Γ from 200 to 500, and the ratio 𝜖𝑒/𝜖𝐵 from 0.1 to 10). In addition, note that 𝑧 has only been measured for a few bursts used in the IC40 analysis, whereas 𝑧≃2 has been assumed for the long bursts with unknown 𝑧. As we illustrated in the previous subsection, this is potentially problematic, which means that a solid lower limit for the prediction can be only obtained for bursts with measured 𝑧 (dashed black curve in Figure 14). Note that our prediction varies not as strongly as one may expect from 𝑓𝜋∝𝑡−1𝑣Γ−4 [20]. The synchrotron losses of the secondaries damp this variation [115]: for smaller Γ, the energy densities in the source increase because of the energy equipartition, and consequently 𝐵′ in (6.6). This increases the energy losses of the secondaries, which means that less energy goes into the neutrinos.

Figure 14 

Numerical prediction (NFC) of the quasi-diffuse flux expected from the bursts used in the IC40-analysis, including the estimates for several model- or method-specific systematical uncertainties (see main text). In addition, the IC40 limit is shown, and two expectations are shown for comparison (IC59+40 from [19] and IC86 extrapolated for 𝐴IC86eff≃3×𝐴IC40eff from IC40; see, e.g., [87]), figure taken from [44].

From the extrapolated IceCube limit for IC86 in Figure 14, it is obvious that IceCube will finally test the parameter space of the GRB fireball models, whereas the current limits already start to enter the meaningful parameter space. On the other hand, it is not clear what values of Γ and 𝑡𝑣 the dominant bursts actually have, and what would contribute most to the neutrino flux. For example, a theoretical study of the dominant Γ to the diffuse neutrino flux has been performed in [30] for different model hypotheses, with a result strongly depending on the model and underlying Γ distribution. In the conventional fireball approach, which is presented above, bursts with Γ lower than 102.5 (as used by IceCube) tend to dominate the expected neutrino flux, such as GRB 091024. In this case, one would expect the prediction rather in the upper half of the shaded region of Figure 14, and a near-future detection in IceCube may be rather likely. On the other hand, recent studies for 𝑡𝑣 seem to prefer larger average 𝑡𝑣≃0.1 s than the 𝑡𝑣≃0.01 s used by IceCube for long bursts [121], which points towards the lower half of the shaded region. Note that in practice, bursts with very large 𝑡𝑣 and Γ (or collision radii) will hardly contribute to the quasi-diffuse limit, which means that the expectation depends on the (weighted) population of bursts with small Γ and 𝑡𝑣 compared to the total population. Independent of these parameters, the neutrino flux prediction is proportional to 1/𝑓𝑒, which means that the final IceCube results will limit the baryonic loading as key parameter.

Finally, note that similar considerations as for the GRBs may apply to AGN models, such as [10, 12, 122], and it is yet to be seen what the impact on the cosmic ray connection is, see [22].

7. Summary and Conclusions

We have discussed the impact of particle physics effects on the neutrino fluxes from cosmic accelerators, based on a generic numerical approach of the photohadronic charged meson production in sources optically thin to neutron escape; compare Figure 6. The starting point has been the proton and photon densities (spectra) within the source, which lead to the secondary meson production. We have included additional 𝑡-channel, higher resonance, and multipion production modes, neutrinos from kaon and neutron decays, the synchrotron cooling of the secondaries (pions, muons, and kaons), the helicity dependence of the muon decays, and the full spectral dependencies. While this approach may not be able to describe every source accurately, it can be regarded as the minimal approach to neutrino production by 𝑝𝛾 interactions using at least the well-known particle physics ingredients, to be extended by further components if applied to particular source classes. We have applied this framework to GRBs in Section 6.1, where the target photons are inferred from the gamma-ray observation. We have demonstrated that it contains all the ingredients in frequently-used analytical models, such as [9, 15, 20, 21], which means that it can be regarded as a numerical refined version of the conventional fireball neutrino model. In addition, we have discussed a generic AGN-like application where the target photons are produced by the synchrotron emission of coaccelerated electrons in Section 5. From the presentation in this paper, it must be clear that this application is at a similar level as the GRB approach in terms of the ingredients. Applications to specific neutrino flux normalization predictions are, however, yet to be performed.

The main consequences of the particle physics effects discussed in this paper can be summarized as follows.(i) Magnetic field effects and flavor mixing change the flavor composition. The assumption of 𝜈𝑒∶𝜈𝜇∶𝜈𝜏 flavor composition at the source of 1 : 2 : 0 (pion beam) can only be justified for 𝐵′≲1 kG, while it will normally change as a function of energy for larger values of 𝐵′.(ii) The neutrino spectral shape for 𝑝𝛾 interactions depends on proton and photon injection spectra and the magnetic field and flavor effects. Roughly speaking, the 𝐸−2𝜈 assumption only holds for the special case (𝜀′)−1 for the photon spectrum. This means that this assumption is too simple for many sources, since the detector response depends on the interplay between spectral shape and differential limit. For example, the detector may miss a neutrino signal because the spectrum does not peak at the right energy, although the fraction of energy going into pion production in the source is the same as in another detectable case.(iii) Additional pion production processes increase the neutrino production significantly and lead to an impact on the spectral shape. Approximations using the Δ-resonance approximation are useful for analytical estimates, but they do not take into account the cross section dependence at high center-of-mass energies. The impact on the neutrino flux is a factor of a few.(iv) Spectral effects, such as the energy dependence of the proton interaction length or the energy dependence of the photon spectrum, leads a significant reduction of the GRB neutrino flux prediction of about one order of magnitude in comparison to analytical estimates. As a consequence, IceCube has just entered the relevant part of the parameter space to test the simplest GRB fireball model. Similar effects on AGN flux models are yet to be tested.(v) Neutrino fluxes from kaon and neutron decays are generically expected in any source. However, the neutrinos from neutron decays show typically up below the peak, the neutrinos from kaon decays above the peak. The relative importance of these effects depends on the strength of 𝐵′ (which may separate the spectra from the different parents) and the maximal proton energy (which controls the neutron decay contribution).(vi) Using the Glashow resonance for the discrimination of 𝑝𝑝 and 𝑝𝛾 interactions in the source is challenging. Especially, it is very hard to infer any physics information from the discovery of Glashow resonant events, whereas the absence of Glashow resonant events (in presence of a neutrino signal) can be interpreted as a 𝑝𝛾 source optically thin to neutron escape.(vii) Already from the current IceCube limits one can extrapolate that conclusions for astrophysical sources will most likely be based on source classes, not individual sources. The aggregation of fluxes from many sources, however, introduces new systematical errors. For instance, the extrapolation from the neutrino fluxes from 117 GRBs to a quasi-diffuse flux implies an error of at least 50% (at the 90% confidence level), from the redshift distribution only.(viii) In the presence of a high statistics neutrino flux (close to the current bound), useful information on possible new physics effects in the neutrino propagation may be inferred from the ratio between muon tracks and cascades. In this case, the energy dependence of the flavor ratio at the source and the energy dependence of the new physics effects may provide the necessary information to identify the effect.

In conclusion, in the presence of data, neutrino astrophysics is now at the point where rough analytical estimates for neutrino fluxes are not sufficient anymore, since in that case no reliable conclusions can be drawn for the astrophysical models and thus the origin of the cosmic rays. Especially the particle physics effects discussed in this paper can be relatively easily taken into account, and they are well known. We have demonstrated with numerous examples that these effects on flux normalization, spectral shape, flavor composition, and neutrino-antineutrino composition cannot be neglected. For example, a correction of the GRB neutrino flux prediction of one order of magnitude has been identified. Therefore, a reliable treatment of the particle physics of the neutrino production should be the state-of-the-art of any neutrino data interpretation in the multimessenger context.

Acknowledgments

The author would like to thank Philipp Baerwald, Mauricio Bustamante, Svenja Hümmer, and Guey-Lin Lin for useful comments, and Svenja Hümmer for providing some figure from her Ph.D. thesis. This publication was funded by the German Research Foundation (DFG) and the University of Würzburg in the funding program Open Access Publishing. It has been supported by Deutsche Forschungsgemeinschaft, Grants WI 2639/3-1 and WI 2639/4-1.