{'fraction_non_alphanumeric': 0.03487256371814093, 'fraction_numerical': 0.01361319340329835, 'mean_word_length': 4.545560359161955, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Describing the Observed Cosmic Neutrinos by Interactions of Nuclei with Matter

IceCube have observed neutrinos which are presumably of extra-galactic origin. Since specific sources have not yet been identified, we discuss what could be learned from the conceptual point of view. We use a simple model for neutrino production from the interactions between nuclei and matter, and we focus on the description of the spectral shape and flavor composition observed by IceCube. Our main parameters are spectral index, maximal energy, magnetic field, and composition of the accelerated nuclei. We show that a cutoff at PeV energies can be achieved by soft enough spectra, a cutoff of the primary energy, or strong enough magnetic fields. These options, however, are difficult to reconcile with the hypothesis that these neutrinos originate from the same sources as the ultra-high energy cosmic rays. We demonstrate that heavier nuclei accelerated in the sources may be a possible way out if the maximal energy scales appropriately with the mass number of the nuclei. In this scenario, neutrino observations can actually be used to test the UHECR acceleration mechanism. We also emphasize the need for a volume upgrade of the IceCube detector for future precision physics, for which the flavor information becomes a statistical meaningful model discriminator as qualitatively new ingredient.

IceCube have observed neutrinos which are presumably of extra-galactic origin. Since specific sources have not yet been identified, we discuss what could be learned from the conceptual point of view. We use a simple model for neutrino production from the interactions between nuclei and matter, and we focus on the description of the spectral shape and flavor composition observed by IceCube. Our main parameters are spectral index, maximal energy, magnetic field, and composition of the accelerated nuclei. We show that a cutoff at PeV energies can be achieved by soft enough spectra, a cutoff of the primary energy, or strong enough magnetic fields. These options, however, are difficult to reconcile with the hypothesis that these neutrinos originate from the same sources as the ultra-high energy cosmic rays. We demonstrate that heavier nuclei accelerated in the sources may be a possible way out if the maximal energy scales appropriately with the mass number of the nuclei. In this scenario, neutrino observations can actually be used to test the UHECR acceleration mechanism. We also emphasize the need for a volume upgrade of the IceCube detector for future precision physics, for which the flavor information becomes a statistical meaningful model discriminator as qualitatively new ingredient.

I. INTRODUCTION
One of the major breakthroughs in astronomy has been the discovery of high energy neutrinos by the IceCube detector at the South Pole. The story started with two PeV neutrino events announced in 2012 [1] as a consequence of changing the search strategy, over 28 events with deposited energies greater than 30 TeV announced in 2013 [2], to 37 events in the current three year data analysis [3]. While about 15 of these 37 events are expected from the backgrounds of atmospheric muons and neutrinos, the excess over the background currently constitutes more than 5σ. Clearly, high-energy neutrino astronomy is emerging as a new discipline with IceCube at the forefront, which is expected to collect 100-200 events (within the current analysis scheme) during its lifetime. While this is significant statistics, it may not be sufficient for resolving individual sources [4] and for precision studies of spectrum and flavor composition. Therefore, the next generation experiments are being discussed, such as a volume upgrade of the IceCube detector (IceCube highenergy extension HEX) and the KM3NeT experiment in the Mediterranean [5]. The results of the IceCube experiment during this and the coming years will have to serve as input for the optimization of these future options.
On the theoretical side, there have been many speculations where these neutrinos would be coming from, see e.g. Ref. [6] for a recent review. It is probably fair to say that there is no general answer on that question yet. Instead, the current state-of-the-art can be recast in a number of conceptual questions: • What is the role of atmospheric neutrinos, especially prompt neutrinos?
• Are there any directional or time-wise clusters, or is the flux isotropic? Are there any correlations with known objects or events?
• Are some of the events of Galactic origin?
• Why are there no events above a few PeV?
• Can the neutrinos stem from the sources of the ultra-high energy cosmic rays?
• Are the observed neutrinos coming from one source class or more? Which ones?
• Is the flavor composition what is expected, or are there deviations indicative for either physics beyond the Standard Model or non-conventional compositions at the source?
• Is there a particle physics origin of these neutrinos, such as dark matter?
In this study, we focus on the conceptual interpretation of the observed events in terms of the spectral shape of the observed flux. That is, we assume that the neutrinos come from one source population with similar properties which is of extra-galactic origin, as there is not yet any evidence for directional clusters. We furthermore assume that the source population is cosmologically distributed such that it roughly follows the star formation rate. We postulate that the neutrinos are produced from interactions between nuclei and matter. This is, in a way, the simplest possible class of models, as the neutrino spectrum directly follows the non-thermal spectrum of the accelerated nuclei; for interactions with radiation, the obtained spectral shape of the neutrinos depends on both the spectra of the interacting nuclei and the target photons, see e.g. Ref. [7] for target photons produced by synchrotron radiation of co-accelerated electrons. However, we do take into account magnetic field effects on the secondary muons and pions, which can significantly alter the neutrino spectra and flavor composition, see e.g. Refs. [8][9][10] especially in the context of gamma-ray bursts [11][12][13] and microquasars [14,15]; see Ref. [16] for a review. While most of these studies discuss interactions between protons and photons, the secondaries produced by interactions between nuclei and matter will be  affected by magnetic fields as well, see e.g. Ref. [14]. In the context of a possible cutoff at PeV energies in recent IceCube observations, magnetic field effects on the secondaries may be a way to decouple the maximal proton from the maximal neutrino energy, see Ref. [7] for a more detailed discussion. Finally, there seems to be condensing evidence for a heavier composition of the UHECRs [17]. We therefore take into account the composition of the nuclei in the sources. We especially discuss if the cutoff at PeV neutrino energies can be consistent with the UHECR paradigm if heavier nuclei are accelerated to higher energies within the sources. In fact, we will demonstrate that one can learn something about the acceleration mechanism in that scenario. Note that the interpretation of the obtained neutrino flux normalization in terms of source luminosity and column depth will be discussed elsewhere.

II. MODEL AND METHODS
The interaction model used in this work is based on the Kelner et al. [18] parameterization for proton-proton interactions, where we take into account the charged pion production explicitly to allow for secondary cooling. For the extension to heavier nuclei Ap interactions, see Ref. [19]. The secondary production Q π [cm −3 s −1 GeV −1 ] is given from the non-thermal nucleon density in the source N A [cm −3 GeV −1 ] and the target nucleon density n p [cm −3 ] by where x A = x/A is the fraction of the nucleus' energy going into the neutrino and f are the scaling functions from Ref. [18] (SIBYLL-based versions) and σ Ap = A 3/4 × σ pp [20]. If the target material is heavier than hydrogen, one may superimpose the nuclei n p A × n A . There are, however, corrections to that, but, as we do not discuss the normalization in this study, that does not affect our results. Note that Eq. (1) can be re-written in terms of the column density L n p , where the interpretation of L in terms of the size of the interaction region depends on the scenario. For example, injecting nuclei into the interaction region with the rate Q A , one can estimate that N A Q A t esc in the absence of disintegration and cooling, and therefore L n p = c t esc n p is the column density which determines the normalization of Eq. (1). The pion and consequent muon decays are computed in the usual way including the helicity dependence of the muon decays, see Ref. [10]. Flavor mixing is taken into account with the best-fit values from Ref. [21] (first octant solution).
The main parameters of the model are listed in Tab. I.
where α is the initial spectral index, which is expected to be α ∼ 2 from Fermi shock acceleration. The maximal energy E max is typically obtained from equating the acceleration timescale with the dominant escape and energy loss timescales in a specific scenario. Since this derivation is highly model-dependent, we keep E max as a model parameter. Note that E max can also be used to simulate a spectral break in the initial spectrum, as it may come from an energy-dependent escape time frequently discussed for starburst galaxies [22], see also Refs. [23][24][25]. In addition, in Eq. (1), N A already corresponds to the result including disintegration and other cooling and escape processes. That means that E max could also describe a spectral break from a cooling process or photodisintegration. The secondary muons and pions are assumed to undergo synchrotron losses and decay governed by the magnetic field B, which impacts spectral shape and flavor composition see Refs. [10,13]. Note that there could be other cooling or escape processes affecting the secondaries, such as adiabatic cooling (see e.g. [26]) or reacceleration [27][28][29]. These effects are, however, modeldependent and typically not the dominant ones shaping the neutrino spectra. Finally, we have the composition A as parameter. In fact, Eq. (1) allows to use an energydependent (average) composition A(E A ) = A(E ν /x A ), which we will use below. Note that for a power law with α = −2, the composition would not affect our results. However, both the cutoff and varying composition will change that conclusion.
For the sake of simplicity, we furthermore assume that the sources do not have large Doppler factors, and that they are cosmologically distributed following the star formation rate by Hopkins and Beacom [30], for steady sources with H(z) the source density normalized to the local source density (H(0) ≡ 1), Ω m = 0.27 and Ω Λ = 0.73. For transients, there will be another factor 1 + z in the denominator of the integrand; we find however that the results depend very little on the details of the cosmological source distribution.
For the fit, we follow Ref. [7], using the up-to-date three year data from Ref. [3]. We use eight bins: four in the reconstructed neutrino energy 30 to 200 TeV, 200 TeV to 1 PeV, 1 to 3 PeV, and 3 PeV to 100 PeV, and one bin for muon tracks and cascades within each energy slot. For the sake of simplicity, we assume that the electromagnetic equivalent energy is roughly 25% of the incident neutrino energy for a muon track, and 75% for a cascade [31]; cascades from neutral current interactions are assumed to be suppressed by the cross sections and their reconstruction at lower energies. The (Poissonian) χ 2 is obtained by bin-wise comparing the observed 36 events, for which energy information is available, with the prediction. The prediction with a free overall normalization is obtained from folding the three flavored neutrino fluxes for one set of model parameters with the corresponding exposures derived from the effective areas in Ref. [2]. Then the atmospheric backgrounds are added. The atmospheric neutrino background is derived from the IceCube observation of muon neutrinos [32,33], and for the atmospheric muons we assume the same shape. The measurement of the atmospheric electron neutrino background is much more uncertain. We therefore extrapolate it from the muon neutrino background, roughly consistent with the flavor composition in [34]. We obtain 3.6 background muon tracks from neutrinos, 3.2 cascades from neutrinos, and 8.6 muon tracks from atmospheric muons, matching the publically available information on the IceCube analysis. This means that the predicted number of muon tracks is slightly higher than the observation (eight), a fact which discussed in detail in Ref. [35]. While this unavoidable tension increases the minimal χ 2 of our fit, it hardly affects the ∆χ 2 , and therefore if of little relevance for the results in this study. As a final step the χ 2 between predicted and observed rates are summed over all bins and minimized over the free normalization of the astrophysical flux to obtain the best-fit for the chosen parameters.

III. RESULTS FOR PROTON-MATTER INTERACTIONS
In Fig. 1, the allowed fit region (filled contours) are shown for in a three parameter (α, B, E max ) model for protons only, where one of the parameters is fixed in each panel. Note that we have chosen a reduced (deposited) energy range 60 TeV ≤ E dep ≤ 3 PeV to compare with the results in Ref. [3]. Most noteworthy, for B ≤ 10 2 G, we obtain a spectral index α = −2.3 ± 0.3 in consistency with the collaboration results. Thus, although our procedure qualitatively deviates from internal analyses of the IceCube collaboration in a few ways (e.g., mapping from deposited to reconstructed energy/energy reconstruction, details on background model, systematical errors), we can roughly reproduce their results in order to test more complicated models. For instance, one can read off from this figure that E max and B can produce a cutoff as alternative to softer spectra. For instance, strong enough B, a spectral index α = −2 is allowed. Note that in the following, we will use the full set of events to make the full use of statistics.
An interesting observation was made by Murase, Ahlers, and Lacki (MAL) [36]: the production of gammarays from π 0 decays in the sources, which are co-produced with the charged pions, may violate the Fermi isotropic background measurements [37]. The highest energy data points are at about 100 GeV, which is significantly below the measured neutrino energies and requires some extrapolation of the spectrum. In addition, gamma-rays may come from higher energies by the initiated electromagnetic cascade. While the details are somewhat modeldependent, we compute the gamma-ray flux injected at the sources based on Ref. [18]. We add a penalty χ 2 , im- , an upper bound. The effect of this penalty can be seen as dashed curves in Fig. 1: it leads to a lower cutoff α −2.2 (1σ). This result is consistent with Ref. [36], where α −2.18 was found. Note that for the chosen star formation evolution of the sources, the main contribution will come from z ∼ 1, where the optical depth at 100 GeV is still small enough such that most gammarays can reach us without being attenuated in photon background fields during their propagation. However, the gamma-ray constraint can be at least partially avoided if a spectral break or lower cutoff in the energy spectrum of the non-thermal nucleons is introduced. We therefore discuss it separately in this study.
Let us now relax the constraint using only the events in the range 60 TeV ≤ E dep ≤ 3 PeV, i.e., we take into account all 36 events with energy information. The fit result is shown in Fig. 2 for a proton composition, where one of the parameters is fixed in each panel (parameter space "section"). The minimal χ 2 is about 9 in the leftpanels, and the χ 2 /d.o.f. is about two. This relatively large value comes from the above mentioned tension between muon track prediction and observation. However, given the small number of bins (eight), it should not be over-emphasized. Comparing to Fig. 1, we note that the spectral index α shifts to softer values (α = −2.7 ± 0.2, 1σ for 1 d.o.f.), see upper left panel for small B. Apart from a better matching at low energies, the information beyond 3 PeV leads to stronger constraints because no neutrinos have been seen there. For the same reason, we have a clear limit in the E max -B plane (upper right panel), where the lower right corner is excluded because it would produce too many high-E events. In that panel, two distinctive regions appear at the 1σ confidence level: one can either produce the cutoff with 10 7 GeV E max 10 8 GeV, or with B ∼ 10 4 G. It is noteworthy that these regions can be potentially discriminated by the flavor composition of the neutrinos: Roughly on the r.h.s. of the dashed-dotted line, the neutrino production will be dominated by pion decays at the highest energies, whereas the muons lose energy faster than they decay ("muon damped source"). As a consequence, only muon neutrinos and antineutrinos will be produced at the source, which leads to a deviation from the canonical (ν e : ν µ : ν τ ) ∼ (1 : 1 : 1) flavor composition at the detector including flavor mixing [38]; see Ref. [16] for a review.
A potential theoretical constraint comes from proton synchrotron losses: not necessarily all of the regions Here only events used by the IceCube collaboration fit [3] in terms of deposited energy 60 TeV ≤ E dep ≤ 3 PeV have been used. Dashed contours refer to including a generic bound by Murase, Ahlers, Lacki (MAL) [36] on the γ-ray emission from π 0 s co-produced with charged pions if the spectrum extends down to 100 GeV, where it has to obey the Fermi isotropic background bound [37]. The best-fit (for the solid contours) is marked by the dot.
shown in Fig. 2 can be reached, as protons may lose energy in magnetic fields faster than they can be accelerated. This can be quantified using an acceleration rate for shock acceleration [39] which corresponds to a constant fractional energy gain per cycle η. For efficient acceleration, one typically assumes η 1. Synchrotron losses are governed by which means that they take over at high enough energies. Equating Eq. (2) with Eq. (3), one obtains where in the latter step m ∝ A ∝ Z was assumed, which is a good approximation for elements heavier than hydrogen. For hydrogen and η = 1, the region where synchrotron losses dominate is on the r.h.s. of the dashed line in Fig. 2, upper right panel. That is, for shock acceleration and and moderately efficient acceleration, that region cannot be reached. We will consider the impact of this theoretical constraint on the fit below. Note that we do not assume relativistic boosting here. For example for gamma-ray bursts, E max 10 9 GeV can be reached for 100 kG in the shock rest frame, which translates into E max ∼ 10 11 GeV in the observer's frame.
So what kind of options to we have to describe the data? In order to illustrate that, three test points in the 1σ region are marked in Fig. 2. We show the obtained best-fit spectra for these test points and electron neutrinos in Fig. 3, upper panel. Furthermore, we show the flavor ratio of muon to electron and tau neutrinos at the detector, which corresponds to the ratio between induced muon tracks and cascades (without efficiencies), in the lower panel. We can identify three options: Point 1. High E max are allowed together with strong magnetic fields. The magnetic fields lead to a cut-off and some characteristic wiggles in the spectrum, which come together with a change of the flavor composition at PeV energies from pion beam to muon damped source, see lower panel of Fig. 3.

Point 2.
For small B, the cutoff can be achieved by an appropriate maximal proton energy, as discussed above.
Point 3. Alternatively, a soft enough spectrum can describe data for small B and large E max ; see also Ref. [40].
We will discuss Point 4 in the next section. Note that in all cases, the normalization (which is a result of the fit) is about 1.5 10 −8 GeV cm −2 s −1 sr −1 at 30 TeV. It is, of course, interesting to discuss what this information tells us about the sources. Point 1 corresponds to sources with strong magnetic fields, such as low luminosity gamma-ray bursts [41], "chocked" gamma-ray bursts [42][43][44], or (extra-galactic) micoquasars or pulsars. Point 2 may correspond to starburst galaxies [22], galaxy clusters/groups [36] or radio galaxies [45]. And Point 3 may come from any extra-galactic population, where the main challenge is to accommodate α −2.2 with the theory of Fermi acceleration. On possibility is the effect of turbulence on Fermi shock acceleration which may cause such effects [46], another one is that the overall spectral index comes from convoluting a harder spectrum with an appropriate luminosity distribution function [47].
Maybe even more interesting is the conceptual question if these neutrinos can come from the sources of the UHECRs. Considering Point 1, which is taking into account the synchrotron loss constraint in Fig. 2, the maximal energy can only be high enough to reach the UHECR range E > 10 10 GeV if large Lorentz boosts are involved. Point 2, on the other hand, cannot be accommodated with the UHECR paradigm, because E max is too low. For Point 3, the discussion is much more complicated: While it can be in principle accommodated with the UHECR paradigm, the soft spectrum tends to lead to neutrino overproduction at PeV energies if one normalizes the UHECR range to normalization. This is discussed for gamma-ray bursts in Refs. [48,49], and, in a more generic context, in Ref. [50]. A model-independent "proof" seems, however, more difficult.
In Fig. 4, we take into account multi-parameter correlations. That is, we show the projections including the minimization of the parameter not shown in each panel. This increases the fit region regions dramatically. On the other hand, we include the Murase-Ahlers-Lacki bound, which leads to α −2.2, and the synchrotron loss constraint, explicitly shown in the upper right panel. These theoretical constraints reduce the size of the fit regions. The main result, which can be read off from the left panels, is that spectral indices compatible with Fermi acceleration are preferred in combination with a cutoff of the maximal proton energy, whereas a wide range of magnetic fields are possible. The best-fit at α 2 and E max 2 10 9 GeV points towards conventional scenarios of Fermi shock acceleration together with a maximal proton energy cutoff compatible with the ankle of the cosmic ray spectrum observed in our Galaxy. It is therefore plausible that the neutrinos are produced under similar conditions, such as in starburst galaxies.

IV. NUCLEI-MATTER INTERACTIONS, AND THE UHECR PARADIGM
Adding the composition to the parameters of the model increases the complexity to a level such that no meaningful information can be obtained from current data due to limited statistics. We therefore focus on the key issue, the proton composition has not been convincingly successful to describe: can the potential cutoff at PeV energies be reconciled with the UHECR paradigm, taking into account that the composition could be as heavy as iron at the highest energies [17]?
Let us assume that the magnetic fields are small enough such that magnetic field effects on the secondaries are negligible. That is a necessary condition such that the maximal energies are not suppressed by synchrotron losses in our standard scenario (apply, for instance, Eq. (4) to dashed curve in Fig. 2 for E max = 10 11 GeV). Let us furthermore assume that α −2 in consistency with the argument in Ref. [50], and that Here the full parameter degeneracy is taken into account, i.e., the χ 2 is minimized over the third parameter in each panel.
Here theoretical exclusion limits are included as χ 2 penalties: the Murase-Ahlers-Lacki bound on γ-ray observations, and the region unreachable because of synchrotron losses dominating the maximal proton energy as a cutoff.
A(E max ) = 56 (iron). We then assume that the maximal energy is element-dependent, which leads to an energydependent composition parameterized as i.e., A(E) ≥ 1 and A(E max ) = 56. The coefficient β describes how A scales with energy, and we treat it as a continuous parameter. It is, however, useful to consider a few examples: Rigidity scaling. This is the most often used approach, also known as necessary condition formulated by Hillas [39]: the Larmor radius has to be smaller than the acceleration region. This can be also described using Eq. (2) by t −1 acc = t −1 lim c/R, where R is the size of the region and t −1 lim the limiting timescale determined by the size of the region (typically the dynamical timescale, escape timescale, or adiabatic cooling timescale). As a consequence, E/Z (the rigidity) is constant for constant B and R, which means that higher energies can be reach for higher charges. Since Z ∼ A/2, one has β 1.
Synchrotron-loss dominated E max . From Eq. (4) describing a shock acceleration scenario, we can immediately read off that A ∝ E 2 , i.e., β = 2, if the maximal energy is limited by synchrotron losses.
It is generically difficult to obtain coefficients β < 1 unless the timescale constraining the maximal energy slightly drops with energy and does not scale with rigidity. This may be achieved in scenarios where e.g. photodisintegration dominates the highest energies [41]. The current best-fit region in terms of β and E max is shown in Fig. 5 for fixed α = −2. It is clear from the figure that E max 10 10 GeV, required to describe UHECR observations, implies that β < 1. Extremely high energies are allowed for 0.05 β 0.35, which requires unconventional assumptions for the acceleration-radiation scenarios -as discussed above. This scenario neither requires strong enough magnetic fields, nor a spectral index softer than α = −2; the cutoff is instead produced by a change of the composition. Note that neutrino data can in that case be used to infer on the acceleration of the heavier elements itself, or to model the injected UHECR composition from the sources in propagation codes. The spectrum corresponding to test point 4 is also shown in Fig. 3 (upper panel). It peaks at somewhat higher energies than the other spectra. It is probably noteworthy that, because of α = −2, there are no issues with the MAL-bound in this case, as we have explicitly tested.

V. FUTURE EXPECTATIONS FOR ICECUBE PERFORMANCE
We finally discuss what can be learned from future upgrades of IceCube. As an example, let us choose the the 3σ-allowed region as a function of E max and B cor- . "IceCube 2014" corresponds to current data (37 events), "IceCube-86 final" to four times the current exposure, and "IceCube-HEX" to a possible high-energy extension with fourty times the current exposure (factor ten larger mass operated over about a decade). The best-fit is marked by a dot.
responding to Fig. 2, upper right panel (but including the MAL bound). The outer (red) curve shows the constraint from current data, the middle green region the expected result from IceCube-86 over about a decade (four times current exposure), the blue region the expected result from a future high-energy extension (HEX) with about ten times the size of IceCube, operated over about a decade. This figure illustrates that while better information will be available in a few years from now with upcoming data, high precision will require an upgrade of IceCube. In this particular case IceCube-86 cannot discriminate between a cutoff from magnetic field effects or a cutoff in the proton spectrum, corresponding to test points 1 and 2. However, the precision in IceCube-HEX will even allow to exploit the transition in the flavor composition expected by magnetic field effects, and discriminate between these regions. This is evident from the event rates. Consider e.g. Point 1 and the muon track and cascade bins between 1 and 2 PeV. For current statistics, the expectations are 0.6 tracks and 0.9 cascades for these bins, whereas for IceCube-HEX, the expectations are 24 muon tracks and 37 cascades. This leads to a statistical relative error of about 1/ √ 24 20% The muon track to shower ratio is about 20% increased in the muon damped case (coinciding with this energy range for the chosen test point) compared to the pion beam case after flavor mixing, see Fig. 3 (lower panel, between 1 and 2 PeV), which means that the statistics between muon tracks and showers becomes meaningful. This is a signif-icant qualitative advance compared to the full statistics IceCube-86 analysis.

VI. SUMMARY AND DISCUSSION
We have studied the interpretation of IceCube data in the production scenario of nuclei-matter interactions, for which the neutrino spectrum follows the non-thermal spectrum of the nuclei. Compared to earlier studies, we have taken into account possible magnetic field effects on the secondary muons and pions and the composition of the accelerated nuclei. We have especially focused on the reproduction of the spectrum, where the flavor composition (cascades versus muon tracks) has been implied as well. We have essentially identified four different options for the initial spectrum of the protons/nuclei to reproduce current data, which can all avoid the overproduction of events beyond a few PeV: 1. An unbroken power law with α ∼ −2 and magnetic fields B ∼ 10 4 G in the source, leading to magnetic field effects on the secondary muons and pions. This may be realized in certain populations of gamma-ray bursts or (extra-galactic) micoquasars or pulsars.

2.
A power law with α ∼ −2 with a break or cutoff about 10 7 to 10 8 GeV (for protons). This option may be favorable for models of starburst galaxies, galaxy clusters/groups, or radio galaxies.
3. An unbroken power law significantly softer than E −2 , where we find α = −2.7 ± 0.2 including all neutrino events in the fit. This option is the simplest possible one, also mentioned in the recent IceCube three-year analysis [3], where we find slightly softer indices because we do not cut as strongly on the energy range. A major drawback has been pointed out by Murase, Ahlers, and Lacki (MAL) [36]: the spectrum consequently exceeds the isotropic gamma-ray background at lower energies, unless a low-energy break is introduced.
4. An unbroken power law with α ∼ −2 and a flat enough change of the composition of the nonthermal spectra in the source from lighter to heavier elements at the highest energies, which can be potentially as high as 10 12 GeV.
While options 1 to 3 have been identified in similar forms in the literature, option 4 is entirely new. Note that regarding options 1 and 2, similar possibilities are obtained for photohadronic interactions producing neutrinos, for which the spectral shape is determined by both the nuclei and target photon spectra; see Ref. [7]. We have been especially interested in which of the above options can be reconciled with the question that the observed neutrinos stem from the sources of the UHE-CRs, see also discussion in Ref. [51]. While option 2) is obviously incompatible with this assumption, option 1) faces the problem that synchrotron losses typically limit the maximal proton energy, and do not allow for high enough maximal energies (unless strong Lorentz boosts are involved, such as in gamma-ray bursts). Option 3), on the other hand, tends to lead to an overproduction of neutrinos at PeV energies exceeding the current Ice-Cube observations when normalized to the UHECR observations, see Ref. [50] for a generic discussion -while it also violates the MAL bound. We have therefore identified option 4) as the most promising one, especially in the light of recent Auger observation pointing towards a heavier composition at the highest energies. It is generically compatible with the UHECR paradigm, as the acceleration energies can be high enough to describe observations, and it is also compatible with the MAL bound.
In fact, we have demonstrated that neutrino data can be used to test the acceleration mechanism in that scenario. The result from neutrino observations may, on the other hand, serve as an input for the UHECR injection at the sources in UHECR propagation models.
Finally, we have pointed out that future precision measurements will require significant volume upgrades, such as an high-energy extension of IceCube or KM3NeT in the Mediterranean. While this statement is most certainly generically true, we have shown that the current IceCube experiment will not be sufficient to exploit the flavor information from muon tracks versus cascades in a statistically meaningful manner to discriminate between pion beam and muon damped sources. The reason is that the flavor transition is expected at high enough energies to describe the cutoff at PeV energies, where the event rates are very low. A significantly larger (about a factor of ten) volume upgrade would, however, allow for a discrimination between options 1) and 2), where the flavor information will be the qualitatively new ingredient. It is therefore important to optimize an upgrade in a way to preserve the flavor or topology identification capability.