A COMBINED MAXIMUM-LIKELIHOOD ANALYSIS OF THE HIGH-ENERGY ASTROPHYSICAL NEUTRINO FLUX MEASURED WITH ICECUBE

Astrophysical neutrinos⁵³ are created in interactions of high-energy cosmic rays with other massive particles or photons (Gaisser et al. 1995). The neutrinos, being electrically neutral and hence unaffected by cosmic magnetic fields, will travel in straight lines from their point of origin to Earth. If the interactions happen close to the acceleration sites of the cosmic rays, they will thus reveal those. Unlike gamma-rays, which are created in the same processes, the neutrinos are also unlikely to be absorbed during their journey (see, e.g., Learned & Mannheim 2000). Because of these properties, neutrinos are ideal messengers to study the sources of high-energy cosmic rays. Astrophysical neutrinos carry information about these sources in their energy spectrum and flavor composition, even if their individual positions on the sky cannot be resolved yet (Hooper et al. 2003; Lipari et al. 2007; Choubey & Rodejohann 2009; Laha et al. 2013). Candidate sources include active galactic nuclei (e.g., Stecker et al. 1991; Mücke et al. 2003), gamma-ray bursts (e.g., Waxman & Bahcall 1997; Guetta et al. 2004), starburst galaxies (e.g., Loeb & Waxman 2006), and galaxy clusters (e.g., Murase et al. 2008), as well as galactic objects like supernova remnants or pulsar wind nebulae (e.g., Bednarek et al. 2005; Kistler & Beacom 2006; Kappes et al. 2007).

To first order, the energy spectrum of astrophysical neutrinos follows that of the cosmic rays at their acceleration sites. If Fermi shock acceleration is the mechanism responsible, a power-law spectrum E^−γ with γ ≃ 2 is expected (Gaisser 1990), although the details depend on the characteristics of the specific sources (see e.g., Becker 2008). Furthermore, the majority of the astrophysical neutrinos are expected to arise from the decay of pions created in cosmic-ray interactions, i.e., $\pi \to \mu +{\nu }_{\mu }$, followed by $\mu \to e+{\nu }_{e}+{\nu }_{\mu }$. The flavor composition resulting from this decay chain is ν_e:ν_μ:ν_τ = 1:2:0. Taking into account long-baseline neutrino oscillations, the flavor composition at Earth is different, approximately 1:1:1 for this scenario (Learned & Pakvasa 1995; Athar et al. 2006). This first-order model of the energy spectrum and flavor composition has often been used as a benchmark scenario in the past.

Second-order corrections to this benchmark model arise from, e.g., muon energy losses (Kashti & Waxman 2005) and muon acceleration (Klein et al. 2013); these effects can alter both the energy spectrum and flavor composition of the astrophysical neutrino flux. While the energy spectrum is difficult to constrain by general arguments, there are two limiting scenarios for the flavor composition at the sources⁵⁴ : muon-damped sources, in which the high-energy neutrino flux is suppressed due to energy loss processes of the muons, and neutron-beam sources, in which the neutrinos are created from neutron rather than pion decays. The flavor compositions at the source for these (idealized) scenarios are 0:1:0 and 1:0:0, respectively (Lipari et al. 2007). Using the neutrino oscillation parameters from Gonzalez-Garcia et al. (2014, inverted hierarchy), the expected flavor transitions for the three source classes discussed here are 1:2:0 $\to$ 0.93:1.05:1.02 (pion-decay), 0:1:0 $\to$ 0.19:0.43:0.38 (muon-damped), and 1:0:0 $\to$ 0.55:0.19:0.26 (neutron-beam).

The implications of deviations from the benchmark scenario have been discussed by several authors. Assuming extragalactic sources, Winter (2013) argues that a spectral index γ ≳ 2.3 is more easily explained in photohadronic scenarios, i.e., by pγ-interactions as the origin of the neutrinos. Similarly, Murase et al. (2013) point out that if a hadronuclear origin (pp-interactions) in extragalactic sources is assumed, measurements of the diffuse extragalactic gamma-ray background (see Ackermann et al. 2015) imply γ ≲ 2.1–2.2 under certain conditions. However, hadronuclear models that take into account the diffuse gamma-ray background as well as recent IceCube data (Aartsen et al. 2015a), which point toward a softer spectrum, have also been proposed, e.g., by Senno et al. (2015). Other authors invoke spectral arguments to propose that the astrophysical neutrino flux is produced by sources within the Milky Way (Neronov & Semikoz 2016; Gaggero et al. 2015). Finally, the implications of flavor ratios different from the benchmark scenario have been discussed by, e.g., Beacom et al. (2003), Lipari et al. (2007), and Vissani et al. (2013). Bustamante et al. (2015) give an overview over the flavor compositions resulting from various standard and non-standard neutrino production and propagation scenarios.

All relevant backgrounds to searches for high-energy astrophysical neutrinos are created in cosmic ray-induced air showers in the atmosphere of the Earth. Of all the particles created in air showers, only muons and neutrinos can reach the IceCube detector. Atmospheric muons constitute, by far, the most abundant background, triggering the detector at a rate of several kHz. Characteristically, they reach the detector from above the horizon and are first detected on the boundary of the instrumented volume. Atmospheric neutrinos are created at a similar rate, but are detected much less frequently due to their small interaction probabilities. They reach the detector from all directions, in particular also from below the horizon. Atmospheric neutrinos that have passed through the Earth are free of muon background, but also hard to distinguish from astrophysical neutrinos. In contrast, atmospheric neutrinos arriving from above the horizon are often accompanied by atmospheric muons (Gaisser et al. 2014), which is never the case for astrophysical neutrinos.

At lower energies, the atmospheric neutrino flux is dominated by so-called conventional atmospheric neutrinos from the decays of kaons and charged pions. Since these particles are likely to interact with air molecules before they decay, the resulting neutrino flux differs from the original cosmic ray flux in its energy and zenith angle dependence: the energy spectrum is steeper (approximately E^−3.7) and the flux is enhanced toward the horizon (Gaisser 1990). An additional contribution to the atmospheric neutrino flux is expected from the decays of heavy, short-lived hadrons containing a charm or bottom quark. At the energies relevant here, these almost always decay before having the chance to interact, giving rise to a flux of prompt atmospheric neutrinos. This flux is predicted to follow that of the cosmic rays more closely, with an energy spectrum of approximately E^−2.7 and an isotropic zenith angle distribution. It is, however, yet to be conclusively observed and predictions for its magnitude vary strongly (Bugaev et al. 1989; Martin et al. 2003; Enberg et al. 2008; Bhattacharya et al. 2015).

The IceCube detector (Achterberg et al. 2006) consists of 5160 optical modules, installed in the ice underneath the geographic South Pole at depths between 1450 and 2450 m. The modules are deployed on 86 strings, where the vertical spacing between two modules on a string is approximately 17 m and the horizontal spacing between two strings is approximately 125 m. Eight of the strings are part of a more densely instrumented region in the center of the detector, DeepCore, which serves as a low-energy extension to IceCube (Abbasi et al. 2012). Each optical module is a spherical glass housing that contains a 10'' photomultiplier (Abbasi et al. 2010), together with digitization electronics (Abbasi et al. 2009). The IceCube Neutrino Observatory also includes a cosmic-ray detector called IceTop, which consists of 81 stations that are installed on the surface above the IceCube detector (Abbasi et al. 2013).

IceCube detects neutrino interactions by measuring the Cherenkov radiation that is induced by the secondary particles created in the interaction. The signature of a neutrino interaction in the detector depends on the flavor of the incoming neutrino as well as on the interaction type. So-called tracks arise from charged-current ν_μ interactions, in which a muon is produced together with a hadronic particle shower at the interaction vertex. At the relevant energies, the muon has a range of several kilometers, leading to an elongated, track-like signature. Due to the long lever arm, the directional reconstruction of these events is very precise, with a typical median resolution of better than 1° (Aartsen et al. 2014a). Since the muon will typically leave the detector, and often also enter it from outside, only a lower limit can be determined for the neutrino energy. On the other hand, charged-current ν_e and ν_τ interactions as well as neutral-current interactions of all flavors give rise to showers.⁵⁵ The electromagnetic and hadronic particle showers produced in these interactions have dimensions that are typically much smaller than the detector spacing, so they appear as point sources of light. As a result, the directional reconstruction of these events is worse than for tracks, typically 15° (Aartsen et al. 2014b). For both event topologies, the energy deposited in the detector, which is a lower bound for the neutrino energy, can be reconstructed with a resolution of about 15% (Aartsen et al. 2014b). For a track event, usually the differential energy loss or the total muon energy is estimated. These quantities are correlated with the energy of the neutrino that produces the muon. However, their absolute values should be interpreted with care (Aartsen et al. 2014b). Because the absolute values are unimportant here, the quantities are quoted in arbitrary units in this article.

The construction of the IceCube detector was finished in 2010 December. However, data had already been taken since 2005, after the installation of the first IceCube string. Data used in this analysis were taken with configurations consisting of 40 strings (2008/2009), 59 strings (2009/2010), 79 strings (2010/2011) and the full 86-string configuration (since 2011).

Searches of this kind select upward-going and horizontal ν_μ-induced tracks. The good angular resolution of these events provides excellent discrimination against downward-going muon events, which can be reduced to a negligible level. Because events with interaction vertices outside the instrumented volume can be selected, this strategy allows the selection of large neutrino samples of high purity. It is, however, restricted by design to directions from which only neutrinos can reach the detector, i.e., the Northern Hemisphere and directions close to the horizon. This implies that atmospheric neutrinos are an irreducible background to searches of this type and can only be distinguished from astrophysical neutrinos on a statistical basis. Because atmospheric neutrinos are predominantly muon neutrinos at all relevant energies, a flux of astrophysical neutrinos will only be visible at very high energies ≳100 TeV. Neutrinos of these energies have a non-vanishing probability of being absorbed during their passage through the Earth (Gandhi et al. 1996); searches for track-like events are therefore most sensitive in the direction of the horizon, where the amount of traversed matter is less than in more vertical directions.

Here we use the neutrino event samples from two different event selections following this approach, performed on data taken in different periods. The first one ("T1") uses data taken in 2009/2010 with the 59-string configuration of IceCube while the second one ("T2") uses two years of data, taken between 2010 and 2012 with the 79-string and the 86-string configurations of IceCube. For both searches, energy and zenith angle information are available for the neutrino sample. Note that we do not use the full event sample provided by the analyses, but restrict ourselves to the high-energy tail of the spectrum. This degrades the ability to constrain the atmospheric neutrino spectrum (and thus, indirectly, the astrophysical spectrum) somewhat, but facilitates the consistent treatment of systematic uncertainties (which are difficult to account for in the high-statistics threshold region of the track-like searches) in the combined analysis presented here.

These searches select shower-like events and are thus sensitive to all neutrino flavors. Lacking the good angular resolution of tracks, searches of this type achieve background suppression by selecting events with a spherical topology, by requiring the interaction vertex to be contained within the instrumented volume of IceCube, and by vetoing events in which light is first detected on the outermost layer of optical modules. This means that the sensitive analysis volume is smaller than the instrumented volume, but allows the selection of neutrinos coming from any direction. However, reliable rejection of muon background events is only possible if there is a sufficiently high probability of detecting the incoming muon. This typically leads to a contamination of residual muon background events at low energies, where muons may produce too little light to be rejected efficiently enough.

Event samples selected by two searches for shower-like events are included here, one performed on data taken in 2008/2009 with the 40-string configuration of IceCube ("S1") and one performed on data taken in 2009/2010 with the 59-string configuration of IceCube ("S2"). Note that the event sample of search S2 has been extended to lower energies compared to what was presented in Schönwald et al. (2013); see Appendix B for further details. Neither analysis originally accounted for the atmospheric self-veto effect, i.e., the possibility to reject atmospheric neutrinos due to atmospheric muons that accompany them (Schönert et al. 2009). For the analysis presented here, the results were corrected for this effect based on a calculation by Gaisser et al. (2014); checks were made to ensure that systematic errors in the description of the effect for these event samples do not affect the results obtained here. The search S1 provides two event samples with slightly different selection criteria, a low-energy sample (with deposited energies between 2 and 100 TeV, "S1a") and a high-energy sample (with deposited energies >100 TeV, "S1b"). Both shower searches provide reconstructed energies for their respective event samples.

The two hybrid searches considered here do not make any requirements on the topology of the selected events, which makes them sensitive to any kind of neutrino interaction and to all directions. The first hybrid analysis ("H1") utilizes the outermost layer of optical modules as a veto against incoming muons, similar to the two shower searches S1 and S2 described above. To ensure a sufficiently high veto probability, it is additionally required that more than 6000 photoelectrons, integrated over all optical modules, be recorded. This value corresponds to a threshold of ∼30 TeV in deposited energy. The second hybrid analysis ("H2") employs a more complex veto algorithm that adapts to the energy of the event and additionally searches for early hits that could be caused by an incoming muon. As a result, the threshold for this event selection is significantly lower, ∼1 TeV in deposited energy. A residual muon background contamination is present in the event samples of both searches, H1 and H2.

Analysis H1 uses data taken in 2010–2013 while analysis H2 uses data taken in 2010–2012, both with the 79-string and 86-string configurations of IceCube. For both event samples, energy and zenith angle information are available. Analysis H2 additionally provides data on the topology of the selected events; those in which more than ten photoelectrons can be associated with an outgoing muon track are classified as tracks, all others as showers.⁵⁶ The atmospheric self-veto effect was implemented in the same way as for the shower searches (see above), based on Gaisser et al. (2014).

Data taken in one period are often subjected to multiple analyses, some of which may have overlapping event samples. Specifically, among the analyses considered here, three (T2, H1, H2) use data taken between 2010 and 2012 and have event samples that overlap. For the combined analysis presented here, however, it is important that the different event samples be statistically independent, i.e., that individual events are not selected more than once. To achieve this, we identify the overlap between the samples and remove the corresponding events from all but one of the samples, in both simulation and experimental data. In particular, we remove events that start inside the instrumented volume of IceCube from the event sample of analysis T2 and we consider only events from the event sample of analysis H2 for which fewer than 6000 photoelectrons were recorded. With these adjustments, all event samples combined in the maximum-likelihood analysis are statistically independent.

The relevant background components to all searches considered here are atmospheric neutrinos and atmospheric muons. We use the same models to describe the neutrino background in all of the different searches. The muon background is more specific and is hence modeled individually for each analysis.

3.1.1. Atmospheric Neutrinos

3.1.2. Atmospheric Muons

So far, the sources of the high-energy neutrino flux measured by IceCube have escaped identification (Aartsen et al. 2014e, 2014f). Here, we therefore test different models for the astrophysical component, based on general assumptions about the neutrino flux. These models are listed in Table 3.

The simplest model is the "single power law" model. We assume that the astrophysical flux arrives isotropically from all directions, that its flavor ratio at Earth is ν_e:ν_μ:ν_τ = 1:1:1, and that it can be described by a simple power law of the form

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\nu }=\phi \cdot {(\displaystyle \frac{E}{100\;\mathrm{TeV}})}^{-\gamma }.\end{eqnarray} \tag{ 1 }$$

Φ_ν denotes the all-flavor neutrino flux, ϕ its value at 100 TeV, and γ the spectral index of the power law.

The "differential" model is based on the same assumptions about isotropy and flavor composition, but it models the astrophysical flux with nine independent basis functions, defined in nine logarithmically spaced energy intervals between 10 TeV and 10 PeV. The normalizations ϕ₁−ϕ₉ of the basis functions are free fit parameters, while the energy spectrum in each interval is assumed to be ∝E⁻². This model is similar to the procedure in Aartsen et al. (2014e, 2015a).

In the "north–south" model, we relax the assumption of isotropy and allow for two independent astrophysical neutrino fluxes, one from the Northern and one from the Southern Hemisphere (separated at decl. δ = 0°), both following a spectrum as defined in Equation (1). While this scenario might lack a good astrophysical motivation, it does allow us to describe the flux separately in two hemispheres that are affected by different detector-related systematic effects. We refrain from testing more complex anisotropic models here because, having been selected for diffuse searches, the event samples used in this analysis currently do not contain unblinded R.A. information. We do, however, note that our simple north–south model could be sensitive to certain anisotropic scenarios such as, e.g., an additional astrophysical component from the inner Galaxy.

Finally, the assumptions about the flavor composition are weakened in the "2-flavor" and "3-flavor" models. In the 2-flavor model, the ν_μ and ν_τ fluxes are assumed to be equal at Earth, while the ν_e flux is allowed to deviate. This relation is, to a good approximation, true for any flavor composition at the source if standard neutrino oscillations transform the neutrino flux during propagation.⁵⁸ In the 3-flavor model, the normalizations of the fluxes of all three neutrino flavors are free parameters, allowing us to test for non-standard oscillation scenarios. In both the 2-flavor and 3-flavor models, the fluxes of the individual neutrino flavors are assumed to have the same energy dependence, i.e., ∝E^−γ, where γ is a free parameter.

In the maximum-likelihood method, the agreement between experimental data and the simulated PDFs is established by means of the test statistic described in Section 3.3.1. Section 3.3.2 lists the systematic uncertainties that we account for in the method. Finally, the calculation of p-values for likelihood ratio and goodness-of-fit tests is explained in Sections 3.3.3 and 3.3.4, respectively.

3.3.1. Test Statistic

We use a binned Poisson-likelihood test statistic to compare the experimental data with the model predictions. The general definition for the test statistic is

$$\begin{eqnarray}&&-2\mathrm{ln}L=-2\displaystyle \sum _{i}\mathrm{ln}(\displaystyle \frac{{e}^{-{\nu }_{i}}{\nu }_{i}^{{n}_{i}}}{{n}_{i}!}),\end{eqnarray} \tag{ 2 }$$

where n_i and ν_i denote the number of observed and predicted events in bin i, respectively.

Systematic effects might change the PDFs that are fit to the data and could thus distort the results of the fit. To avoid this, we have parameterized the impact on the PDFs of all relevant systematic effects and included them in the fit procedure as nuisance parameters. For each nuisance parameter θ there is a prior, an additional, Gaussian-shaped penalty term in the likelihood function that penalizes deviations from the default central value θ^*. The width σ[θ] of the prior is based on the uncertainty associated with the systematic effect. With m nuisance parameters included, Equation (2) now reads

$$\begin{eqnarray}&&-2\mathrm{ln}L=-2\displaystyle \sum _{i}\mathrm{ln}(\displaystyle \frac{{e}^{-{\nu }_{i}}{\nu }_{i}^{{n}_{i}}}{{n}_{i}!})+\;\displaystyle \sum _{j=1}^{m}{(\displaystyle \frac{{\theta }_{j}-{\theta }_{j}^{*}}{\sigma [{\theta }_{j}]})}^{2}.\end{eqnarray} \tag{ 3 }$$

The individual systematic effects considered as nuisance parameters are discussed in greater detail in the next section.

3.3.2. Systematic Uncertainties

3.3.3. Likelihood Ratio Tests

Two different models H0 and H1 can be compared with likelihood ratio tests. For this, we use the quantity

$$\begin{eqnarray}&&-2{\rm{\Delta }}\mathrm{ln}L=-2\mathrm{ln}{L}_{{\rm{H}}0}/{L}_{{\rm{H}}1},\end{eqnarray} \tag{ 4 }$$

which Wilk's theorem predicts to be ${\chi }^{2}$-distributed with k degrees of freedom, where k is the difference in the number of parameters between H1 and H0 (Wilks 1938; Olive et al. 2014, chapter 38).⁵⁹ The distribution can deviate from a χ²-distribution if the sample size is small or a parameter is close to a physical bound. In this case, the exact distribution of $-2{\rm{\Delta }}\mathrm{ln}L$ can be computed from toy Monte-Carlo experiments, generated from the best-fit parameter values of model H0. The p-value of the likelihood ratio test is then given by the percentage of toy experiments that have a larger value of $-2{\rm{\Delta }}\mathrm{ln}L$ than observed in data. The p-values of the likelihood ratio test quoted in this article were obtained with this procedure.

We also employ likelihood ratio tests to determine confidence intervals on the model parameters. For each parameter value, we perform a likelihood ratio test between the model with the parameter constrained to this value and the model with the parameter unconstrained. This procedure is known as a profile likelihood scan. In the χ²-approximation that we use here, the 68% and 90% confidence level intervals are given by the values at which $-2{\rm{\Delta }}\mathrm{ln}L=1$ and $-2{\rm{\Delta }}\mathrm{ln}L=2.71$, respectively. Similarly, two-dimensional confidence contours can be obtained by constraining two parameters simultaneously.

3.3.4. Goodness-of-fit Tests

With a slight modification, the quality of the fit of a single model can be assessed using the test statistic. Similar to Equation (4), we define

$$\begin{eqnarray}&&-2\mathrm{ln}L/{L}_{\mathrm{sat}}=-2\displaystyle \sum _{i}\mathrm{ln}({e}^{{n}_{i}-{\nu }_{i}}{({\nu }_{i}/{n}_{i})}^{{n}_{i}}),\end{eqnarray} \tag{ 5 }$$

where L_sat is the likelihood of the model that exactly predicts the observed outcome (i.e., ν_i ≡ n_i). After adding the same term as in Equation (3), this quantity is minimized by the same parameter values that minimize the test statistic defined there. Moreover, in the large-sample limit, the minimum value follows a ${\chi }^{2}$-distribution and as such allows the calculation of a goodness-of-fit p-value (Olive et al. 2014, chapter 38). Since our samples are not large, we determine the distribution of the modified test statistic by generating toy Monte-Carlo experiments from the best-fit parameter values. Comparing this distribution with the value observed in experimental data gives the goodness-of-fit p-value.

The best-fit parameter values for the single power law model are listed in Table 5, including all nuisance parameters. Specifically, the all-flavor astrophysical neutrino flux is

$$\begin{eqnarray}&&\phi =({6.7}_{-1.2}^{+1.1})\cdot {10}^{-18}\;{\mathrm{GeV}}^{-1}\;{{\rm{s}}}^{-1}\;{\mathrm{sr}}^{-1}\;{\mathrm{cm}}^{-2}\end{eqnarray} \tag{ 6 }$$

at 100 TeV and the best-fit power law has a spectral index of

$$\begin{eqnarray}&&\gamma =2.50\pm 0.09.\end{eqnarray} \tag{ 7 }$$

This measurement is valid for neutrino energies between 25 TeV and 2.8 PeV. This energy range was determined by successively removing events, ordered in energy, from the simulation data and repeating the fit with ϕ and γ constrained to their best-fit values; its bounds denote the energies at which the test statistic worsens by $-2{\rm{\Delta }}\mathrm{ln}L=1$.

We obtain a reasonable p-value of 37.6% from a goodness-of-fit test for this model. A power law with a fixed index of γ = 2 is disfavored with a significance of 3.8σ (p = 0.0066%) in a likelihood ratio test with respect to the model with a free spectral index. We also tested a single power law model with a high-energy exponential cut-off as well as a model that consists of two isotropic astrophysical components, each described by a power law. Neither model gave a better description of our data. A power law with a fixed index of γ = 2 and a high-energy exponential cut-off is still disfavored with a significance of 2.1σ (p = 1.7%) with respect to the model with a free spectral index. The best-fit cut-off energy is $({1.6}_{-0.7}^{+1.5})\mathrm{PeV}$ for a fixed spectral index of γ = 2, with a corresponding normalization of $\phi =({5.2}_{-1.5}^{+1.9})\times {10}^{-18}\;{\mathrm{GeV}}^{-1}\;{{\rm{s}}}^{-1}\;{\mathrm{sr}}^{-1}\;{\mathrm{cm}}^{-2}$. No cut-off is fitted for the best-fit spectral index of γ = 2.5.

The likelihood analysis favors a prompt atmospheric component equal to zero, with a 90% C.L. upper limit of 2.1 times the (modified) model prediction by Enberg et al. (2008). This limit is slightly higher than that obtained in Aartsen et al. (2015a); this is due partly to the way the data samples are combined here (cf. Section 2.4) and to the different treatment of the energy scale uncertainty.

The correlation between the parameters ϕ and γ is visualized in Figure 1. The behavior of these parameters as a function of the magnitude of the prompt component is shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the results for the parameters ϕ_prompt, ϕ, and γ again (row "combined"), together with results obtained from the application of the maximum-likelihood method described in this paper to the individual event samples. For the individual fits, the event samples that were reduced to remove overlap (see Section 2.4) were restored to their original size. Furthermore, the normalization of the conventional atmospheric neutrino flux was treated as a nuisance parameter with a prior of 25% around the model prediction for the fits on the event samples of analyses S1, S2, and H1, since this component is not well constrained by these samples alone. A large prompt component instead of an astrophysical component is preferred in the fit of analysis S2 because these two components are close to degenerate if the only observable is the deposited energy and the astrophysical spectrum is steep. Note that some of the results obtained with the individual samples differ slightly, although not significantly, from the originally published results (see references in Table 1). For instance, a softer spectral index than presented in Aartsen et al. (2015c) is obtained for the sample T2 because only the high-energy data are used here; this difference is well within the uncertainty on this parameter. The somewhat harder spectral index and lower normalization for sample H1 measured in Aartsen et al. (2014e) are a result of the reduced energy range (>60 TeV) used in the spectral fit there.

Zoom In Zoom Out Reset image size

Experimental and simulated data, weighted to the best-fit result, are shown in Figure 4 for all event samples included in this analysis. The data are projected onto observables that are used in the likelihood fit. For the hybrid analyses H1 and H2, we show the energy distribution in two different zenith angle bins.

Zoom In Zoom Out Reset image size

Finally, the best-fit spectra for atmospheric and astrophysical neutrinos in the single power law model are shown in Figure 5.

Zoom In Zoom Out Reset image size

Table 6 lists the best-fit parameter values for the differential model. ϕ₁ through ϕ₉ are the normalizations of the corresponding basis functions, each assumed to follow a spectrum ∝E⁻² and defined in logarithmically spaced energy intervals between 10 TeV and 10 PeV. The resulting differential spectrum is shown in Figure 6, together with the single power law model (cf. previous section).

Zoom In Zoom Out Reset image size

The best-fit values for the parameters of the astrophysical neutrino flux, separated into Northern and Southern Hemisphere, are listed in Table 7. Most notably, the best-fit spectral index in the northern sky is ${\gamma }_{{\rm{N}}}={2.0}_{-0.4}^{+0.3}$, while in the southern sky it is ${\gamma }_{{\rm{S}}}=2.56\pm 0.12$. This discrepancy with respect to the single power law model corresponds to a statistical p-value of 13% (1.1σ).

For the 2-flavor model, the best-fit values for the astrophysical ν_e and ν_μ + ν_τ fluxes are listed in Table 8. The constraints on the astrophysical spectral index, which is assumed to be the same for both flavor components, are identical to those obtained in the single power law model. Furthermore, the results on the flavor composition do not depend strongly on the value of this parameter. We measure an electron-neutrino fraction of the astrophysical neutrino flux of 0.18 ± 0.11 at Earth. Figure 7 compares this value to fractions expected for different composition scenarios at the sources of the astrophysical flux.

Zoom In Zoom Out Reset image size

Finally, Table 9 gives the best-fit values for the astrophysical fluxes in the 3-flavor model. As in the 2-flavor model, we find that the results do not change significantly when varying the spectral index of the astrophysical flux within its uncertainties. Note that none of the analyses considered here positively identifies charged-current tau-neutrino interactions, which, depending on the tau decay mode, have a signature very similar to that of charged-current electron-neutrino (∼83%) or muon-neutrino (∼17%) interactions. For this reason, the ν_τ flux is easily conflated with the ν_e and ν_μ components in the 3-flavor model. Figure 8 shows a profile likelihood scan of the flavor composition as measured on Earth, again compared to ratios expected for different source composition scenarios. Performing likelihood ratio tests between these points and the best-fit hypothesis, we find p-values of 55% (0:1:0), 27% (1:2:0), and 0.014% (1:0:0), respectively. Hence, our data are compatible with a pure muon-neutrino composition and the generic pion-decay composition at the source, but reject a composition consisting purely of electron neutrinos at the source with a significance of 3.6σ.

Zoom In Zoom Out Reset image size