# How bright can the brightest neutrino source be?

###### Abstract

After the discovery of extraterrestrial high-energy neutrinos, the next major goal of neutrino telescopes will be identifying astrophysical objects that produce them. The flux of the brightest source , however, cannot be probed by studying the diffuse neutrino intensity. We aim at constraining by adopting a broken power-law flux distribution, a hypothesis supported by observed properties of any generic astrophysical sources. The first estimate of comes from the fact that we can only observe one universe, and hence, the expected number of sources above cannot be too small compared with one. For abundant source classes such as starburst galaxies, this one-source constraint yields a value of that is an order of magnitude lower than the current upper limits from point-source searches. Then we derive upper limits on assuming that the angular power spectrum is consistent with neutrino shot noise yet. We find that the limits obtained with upgoing muon neutrinos in IceCube can already be quite competitive, especially for rare but bright source populations such as blazars. The limits will improve nearly quadratically with exposure, and therefore be even more powerful for the next generation of neutrino telescopes.

###### pacs:

95.85.Ry, 98.70.Vc## I Introduction

IceCube firmly detected astrophysical neutrinos, but currently, it is not possible to identify a neutrino source and the distribution of neutrino events is consistent with being isotropic IceCubeScience ; Aartsen:2014gkd ; Aartsen:2015knd ; IceCubePS ; IceCubeAnis ; Aartsen:2016xlq . Accumulating more and more data of the diffuse intensity will sharpen constraints on an average source flux, but the flux of the brightest source cannot be probed directly with this approach as long as the distribution remains consistent with isotropic. How bright can the brightest neutrino source be? This is the next question that needs to be addressed. Searches for point-like sources determined that the upper limit (post-trial and per neutrino flavor) on the flux of the brightest neutrino source, , ranges from to TeV cm s, depending on declination and assuming energy spectrum IceCubePS .

Here, we address the same question by taking a different approach. In particular, we implement a statistical distribution for the flux of neutrino sources, a more realistic hypothesis than the single-flux population assumed in, e.g., Refs. IceCubePS ; IceCubeAnis . By constraining the shape of the source flux distribution with observables such as the intensity and anisotropies of the diffuse neutrinos, we will derive constraints on .

Our approach is twofold. First, we discuss estimates on that are intrinsic to the fact that we only have access to one universe to sample the source distribution. If the expected number of sources at becomes much smaller than one, then it is unlikely that one could observe larger fluxes in this universe. We show that if the number of sources producing the diffuse neutrino flux measured by IceCube is greater than 10, then this one-source limit of is smaller than the upper limits from Ref. IceCubePS . Thus, our findings allow us to make statements for a flux regime that is still unprobed by IceCube.

Recent analysis of the angular power spectrum found no significant clustering of multiple events IceCubeAnis . As our second approach, we set upper limits on based on this null result, and show that they are tighter than what is inferred from the search for point-like sources, at least for rare source populations. These constraints on are effective in a regime where the one-source limit is above the point-like source limit, showing that the two strategies followed are complementary. We find that the method is particularly constraining even with the current IceCube exposure if we adopt upgoing muon neutrino events Aartsen:2016xlq , which would provide a critical test for blazar interpretation as the origin of the diffuse neutrino flux. We also find that the limits obtained from the angular power spectrum improve quadratically with the exposure. Thus, they provide an extremely powerful probe for the next generation of neutrino telescopes, such as IceCube-Gen2 IceCubeGen2 and KM3NeT KM3NeT .

In this paper, we constrain the flux of the brightest source (rather than, e.g., its joint luminosity and distance), as it is the quantity that is directly relevant to detectability of the neutrino sources—a goal yet to be achieved. Although the flux is a phenomenological quantity, this way, we can make our discussions model independent. Another complementary approach would be to use typical luminosity and density of each source. Although these are more physical quantities, the discussions tend to be highly model dependent. We provide useful conversion formulae for a representative case.

This paper is organized as follows. After introducing relevant formulation of the flux distribution and its relation to the intensity and angular power spectrum in Sec. II, we discuss current constraints on using the one-source argument and the angular power spectrum in Secs. III and IV, respectively. In Sec. V, we apply these generic discussions to several cases of known source populations. Section VI is then devoted to what is expected in the future, before briefly concluding in Sec. VII.

## Ii Formulation

We define as the total number of sources from all sky and as their surface number density. The source flux distribution function is defined as and we also use the equivalent probability density function of the single source . Our hypotheses on the form of are rather mild: We assume that the distribution follows a broken power-law with physically motivated parameters. In particular, denotes the slope of the distribution, , above a characteristic flux . We assume , which is compatible with what is observed in sources detected in other wavelengths such as gamma rays, e.g., blazars BLLacLF ; FSRQLF ; Ajello:2015mfa ; Inoue:2015kuy , star-forming galaxies TAM ; Feyereisen2016 , and radio galaxies DiMauro:2013xta ; Hooper2016 . In fact, if these sources are distributed homogeneously in a local volume where cosmological effects can be ignored (), it is well known that the flux distribution reduces to the Euclidean limit, i.e., Peacock . This is expected, in particular, for the brightest sources (since these are likely to be nearer to us than the fainter members of their source class), and therefore, will be our reference value. For fluxes smaller than , the slope of the distribution must flatten in order to avoid divergences (cf. Olbers’ paradox). We assume for with . The flattening of the slope at low fluxes is, again, supported observationally BLLacLF ; FSRQLF ; Ajello:2015mfa . The top panel of Fig. 1 schematically shows this distribution. A discussion of flux distributions with the assumption on the power-law slopes is postponed until Appendix D.

In a pixel with a size that roughly corresponds to
the angular resolution of the detector, there are on average sources, with .
Then, the flux per pixel is given by the sum of the fluxes of individual sources.^{1}^{1}1In general, is non-integer, and thus a more precise expression is given by a convolution
with a Poisson distribution.
The mean and variance of the flux distribution per pixel, ,
is simply given by times the mean and variance of the
flux distribution per source, :

(1) | |||||

(2) |

where and indicate averages taken over and , respectively. Under our assumptions for , it is straightforward to show that

(3) | |||||

(4) | |||||

where and are both constants of order unity. Note that, in Eq. (4), instead of integrating up to infinity, we truncated at . We define as the typical number of sources per pixel around flux , i.e., , and similarly, we define and corresponding to and , respectively. Then, we obtain the following for the first two moments of the flux distribution:

(5) | |||||

(6) |

Equivalently, the intensity of the neutrino flux (also often referred to as ) and its Poisson angular power spectrum are, respectively,

(7) | |||||

(8) |

The middle and bottom panels of Fig. 1 show the flux distribution multiplied by appropriate powers of such that the area below the curves is proportional to and of , respectively.

In the following, expressions with an explicit index , such as and , represent differential quantities with respect to energy, and those without the index are the quantities integrated over the energy.

## Iii One-source constraint

We are limited to observe a single universe, which then limits our capability to constrain physical quantities. Specifically, we cannot probe arbitrarily large fluxes, because once the number of sources expected at such fluxes becomes smaller than one, it is unlikely to reconstruct the distribution in the region. We define the one-source limit on the flux of the brightest neutrino source, , such that only with a small probability could we find at least one source brighter than in the entire sky.

The mean number of sources above is given by , where is the complementary cumulative distribution function corresponding to . Using the Poisson distribution with this mean, the probability of finding no source brighter than is . By solving this for a power-law , we obtain

(9) |

which further translates into

(10) |

In Eq. (9), depends only on the properties of the source distribution function. In Eq. (10), on the other hand, it is recast in terms of the measured intensity and the free parameter . For the Euclidean case (), . We assume that the intensity refers to neutrinos per flavor, and where necessary, that flavor democracy holds, i.e., . For an assumed energy spectrum (in order to allow a direct comparison with earlier results IceCubePS ), TeV cm s sr, even though a softer spectrum provides a better fit IceCubeICRC .

Figure 2 shows the one-source limits on the flux of the brightest source, , as a function of obtained with
Eq. (10) for a few values of and
. For ease of comparison with the existing literature, these upper limits are presented at 90% confidence
level (CL; ).^{2}^{2}2Taylor expanding for small ,
the reader may approximately rescale these upper limits from a significance to any desired significance with the ratio
.
The upper limit clearly gets weaker when .
For and , Eq. (10)
yields .
For comparison, we also show from Eq. (7)
with its uncertainty from the estimated error on (orange band),
and the upper limit from the search for point-like
sources IceCubePS (horizontal dashed line).
For derivation of the latter, see Appendix A; see also
Ref. Ahlers:2014ioa for an estimate of the sensitivity when the
source density is modeled to follow the star-formation rate.

For source numbers greater than around 10, the one-source limits reach below the upper limit from the search for point-like sources IceCubePS . In other words, finding a source at the flux level close to the point-source upper limits for a source population characterized with (and and ) is unlikely with a chance probability of .

The flux cutoff is caused by either an intrinsic cutoff of the luminosity function or by the volume effect, the latter of which is the case for Euclidean sources (; see Appendix B). Then, Eq. (10) can be regarded as a prediction of . For a given , has to be located between the values of Eq. (10) evaluated with and , at 90% CL. This is shown as a blue band in Fig. 2 for .

We note that it is possible for the modeled population of sources to give only a subdominant contribution to the diffuse neutrino intensity. Indeed, Refs. Bechtol:2015uqb ; Glusenkamp:2015jca ; Murase:2015xka suggest that neither starbursts nor blazars can explain the entirety of the observed neutrino flux. In that case, the one-source constraints become even tighter, as in Eq. (10) should be replaced by , where is the fraction of the measured intensity explained by the source class under investigation. Having in Eq. (10) will improve these limits considerably.

## Iv Angular power spectrum

The maximum flux can also be constrained by measuring the variance of the source flux distribution; this information is essentially equivalent to the angular power spectrum. Indeed, if is too large, only a few of the brightest sources would be enough to make the distribution of neutrinos highly anisotropic by yielding clustered events, in conflict with what is measured IceCubeAnis .

### iv.1 Formalism

The number of neutrino counts per pixel is obtained by multiplying the flux per pixel by the exposure, i.e., the product of the effective area and the live time of the telescope. Note that since the energy spectra of the astrophysical and atmospheric neutrinos differ, so do the corresponding exposures for each component, denoted by and , respectively. The probability distribution of the number of neutrinos per pixel is therefore obtained by convolving the per-pixel flux distribution and the Poisson distribution with mean :

(11) |

where is the flux of the atmospheric backgrounds, which are assumed to be isotropic. It is straightforward to obtain the moments of the distribution of :

(12) | |||||

(13) |

The first term of Eq. (13) corresponds to the Poisson angular power spectrum that originates from discreteness of the sources [Eq. (8)], and the second corresponds to the shot-noise of the neutrinos,

(14) |

where is the surface density of atmospheric background events (see, e.g., Refs AK ; AKNT ; Fornasa2016 in the case of gamma rays).

The rms error for the angular power spectrum at multipole is

(15) |

where is a fractional sky coverage and is a beam window function corresponding to the angular resolution of IceCube AK ; AKNT ; Fornasa2016 . Since the purpose of this study is to obtain a simple estimate of the current limits and future sensitivity rather than accurate values, we assume . Given the null results from the anisotropy analysis IceCubeAnis , we estimate the upper limits on the Poisson angular power spectrum with

(16) |

where (1.64) corresponds to the limits at 90% (95%) CL. By solving this as an equality for , we obtain such that . Then, by using Eqs. (7) and (8), we obtain the corresponding upper limits on as

(17) |

To summarize, our estimates of will rely on observable inputs (, ), instrumental inputs , and theoretical inputs , which we will discuss for different source populations in Sec. V. We present this analysis applied to two of the “clean” datasets of high-energy neutrinos from IceCube.

### iv.2 High-Energy Starting Events (HESE)

Since we care about the angular power spectrum of astrophysical sources, we consider in the first instance only the High-Energy Starting Events (HESE) dataset IceCubeICRC , a relatively clean event sample consisting of showers and contained tracks at the highest energies.

We estimate by using (39) and four years of IceCube exposure for the muon (electron and tau) neutrinos for the tracks (showers), a full-sky coverage , the energy-dependent HESE effective area (from 1 TeV to 10 PeV) from Ref. IceCubeScience , and the live time of the telescope (taken accordingly to be 1347 days). The expected number of neutrinos is consistent with the results of the four-year searches from Ref. IceCubeICRC : For an energy spectrum proportional to , we find the total number of neutrinos . The rest of the measured events should be attributed to atmospheric backgrounds and statistical fluctuations. We also adopt angular resolutions of the order of the median angular resolution of the HESE events, namely and 20 for tracks and showers respectively.

With these parameters, we obtain an upper limit on the Poisson angular power spectrum of

(18) |

for the HESE tracks and

(19) |

for the HESE showers. Since the track events provide tighter constraints by more than one order of magnitude, in the following, we will focus only on the flux limits due to the tracks, so the intensity used in Eq. (17) is that of the muon flavor.

Figure 3 shows the derived from HESE tracks, as a function of and for different values of and . Values of larger than the solid or dotted lines are excluded, as the term due to the flux variance in Eq. (13) would have been detected in Ref. IceCubeAnis . For small values of (at most below 50, in the case with and ), the upper limits obtained here are more stringent than those by the search for point-like sources IceCubePS , let alone the one-source constraints considered earlier. Note, however, that this upper limit is based on the assumption that ; otherwise the source flux distribution would be proportional to with a truncation at (see Appendix D).

### iv.3 Upgoing muon neutrinos

It is possible to repeat the analysis above for high-energy upgoing tracks, for which rather than requiring the interaction vertex be contained one uses the Earth itself as a veto against atmospheric muon backgrounds Aartsen:2016xlq . Above 300 TeV, it is possible to estimate using the best-fit powerlaw models of astrophysical flux Aartsen:2016xlq and the conventional atmospheric background Honda:2015fha . We adopt a sky coverage of , as well as the energy-dependent effective area and construction-dependent livetimes of the telescope from Ref. Aartsen:2016xlq . This corresponds to and , and is consistent with Fig. 1 from Ref. Aartsen:2016xlq where a cursory inspection yields roughly 60 and 10 events above 300 TeV respectively. We adopt an angular resolution of , better than for the contained events of the previous section since the outermost optical modules of IceCube are used to improve pointing rather than as a veto. With these parameters, we obtain an upper limit on the Poisson angular power spectrum of

(20) |

from uncontained, upgoing tracks above 300 TeV.

Figure 4 shows the derived from upgoing tracks. These limits are many orders of magnitude stronger than the limits from HESE as a result of the improved angular resolution and the much larger exposure. The “pivot point” for which the limit is independent of is also below the point-source searches. In addition to these upper limits, we show the region containing the brightest sources at 90% CL derived in Sec. III. The absence of anisotropies will clearly constrain rare sources better than point-source searches for . Complementarily, for more abundant sources, the point-sources searches do not cut into the brightest-source containment band, so we should not expect (with 90% CL) to have seen them yet anyway. This is especially true if we expect multiple source populations to contribute to this flux, since for populations contributing fractions of the isotropic flux this band is even lower. Even allowing for uncertainties in , these two complemetary constraints (which rely only on the physically-motivated assumption that source fluxes are power-law distributed) jointly place a stronger constraint on the brightness of the brightest high-energy neutrino source than current point-source searches.

## V Application to known source populations

Although we aim to make our discussion as generic as possible, such that it can be applied even to unknown classes of astrophysical sources that may contribute at high energies Ando:2017alx , it is certainly of interest to discuss known source populations in this context. We discuss mainly two source classes commonly thought to be the origin of the observed isotropic flux: BL Lacs Giacinti:2015pya ; Righi:2016kio ; Padovani:2015mba ; Padovani:2016wwn and starburst galaxies Murase:2013rfa ; TAM ; Anchordoqui:2014yva ; Chang:2014hua ; Senno:2015tra ; ATZ .

### v.1 Phenomenological representation

The phenomenological parameterisation of a source population we introduced in Sec. II can be summarised by the tuple . The parameters for sources from the second catalog of hard Fermi sources (2FHL; mostly BL Lacs) and starburst galaxies are (2.5, 1.7, ) 2FHL and (2.5, 1.0, ) Feyereisen2016 , respectively. These are estimated from their gamma-ray observations (with help of infrared observations in the case of the starbursts) and assuming a linear correlation between the gamma-ray and neutrino luminosities, . This is well supported for the case of starbursts, which emit neutrinos through interaction TAM ; ATZ . For the blazars emitting through interaction, on the other hand, the relation between the gamma-ray and neutrino luminosities is more complicated and model dependent, but see, e.g., Ref. Padovani:2015mba for a model of linear scaling. Other cases with stronger dependence can also be accommodated with similar parameters: e.g., for the BL Lacs with scaling Tavecchio2015 , and for the flat-spectrum radio quasars with MID . See Appendix C for more discussions for these cases.

With these parameters, Figs. 2 and 4 show that the 90% CL upper limits on the flux of the brightest high-energy neutrino source, are

(21) |

for the 2FHL sources, based on the angular power spectrum constraint, and

(22) |

for the starbursts, based on the one-source constraint. Recall that these upper limits are on the flux per flavor of a population contributing a fraction of the observed astrophysical flux, assuming an spectrum, and requiring (for the former constraint) an absence of detectable anisotropies.

### v.2 Physical representation

Up to this point, we considered , and as free parameters. Another complementary representation is to use more physical quantities such as luminosity and density of the sources, although the discussion will be model dependent. The latter approach was taken in, e.g., Refs. Silvestri:2009xb ; Ahlers:2014ioa ; MuraseWaxman , where sources were assumed to have the same luminosity. These two representations can be converted from one to the other through

(23) | |||||

(24) | |||||

(25) | |||||

in the case of . Typically and for the starbursts and BL Lacs, respectively MuraseWaxman . However, these relations apply only to mono-luminous case as was studied in the literature. See Appendix B for their derivation and more discussions.

In Fig. 4 (and those that follow), we show reference fluxes of some well known sources for each class: Mkn 421 for the BL Lac blazars and M82 or NGC 253 for the starbursts. Mkn 421 is predicted to have a flux around in a model of Ref. Tavecchio2015 . For M82 and NGC 253, we estimate the neutrino luminosity from the gamma-ray luminosity of these sources M82NGC253 , and then by converting to the neutrino luminosity assuming interaction TAM . In addition, we show predicted neutrino flux from the most promising radio galaxy, Cen A, assuming production from interaction Hooper2016b . We assume that these sources are drawn from a population of emitters with the same luminosity. Thus, the number of sources can be estimated by Eq. (24) with , , and typical neutrino luminosity for this population found in Ref. MuraseWaxman .

### v.3 Discussion

All these sources fall within the 90% region of predicted with the one-source argument with (shown as a red band in Fig. 4) and so a source from any of these populations is plausibly the brightest neutrino source. A slight tension exists for Mkn 421, but Ref. Tavecchio2015 predicts several more BL Lacs with similar flux such as PKS 2155-304, and the tension might go away when using a fraction for the blazars. The 90% containment band for is an order of magntidue below the point-source constraint, suggesting it would be unlikely to identify starburst galaxies amongst the brightest neutrino sources. This result is consistent with the analyses in Refs. Feyereisen2016 ; MuraseWaxman .

The angular power spectrum is especially constraining for rare sources such as blazars. The upper limit, Eq. (21), is nearly an order of magnitude lower than the 90% containment band for and the predicted neutrino flux of Mkn 421. The isotropy of the upgoing flux, if confirmed with the current IceCube exposure, will force us to abandon the assumption that they contribute a fraction of the high-energy neutrino flux. This not only eases the aforementioned one-source tension for Mkn 421, but furthermore is consistent with the analysis in Ref. Glusenkamp:2015jca .

## Vi Prospects for the future

In this section, we forecast the prospects for studying the flux of the brightest source with the next generation of neutrino telescopes, under the assumption that anisotropy searches will continue to yield null results in the future.

The angular power spectrum will become much more powerful for IceCube-Gen2 IceCubeGen2 and KM3NeT KM3NeT . This is because of the strong dependence of on from Eq. (17), where improves with exposure as described by Eq. (15). For Euclidean sources (), the upper limit improves quadratically with exposure: . The anticipated tenfold increase in exposure expected for IceCube-Gen2 with respect to the current IceCube IceCubeGen2 will yield hundredfold improvement on if the observed angular power spectrum remains consistent with isotropy, before even accounting for any improvements in angular resolution.

Figures 5 and 6 summarize future prospects for upper limits on the flux of the brightest source, drawn from a population described by and , with an improved track angular resolution and larger exposures than acheived today (cf. Table 1). For comparison, we scale down the upper limit from the search of point-like sources by a factor of , assuming that these analyses are already background limited; the value of from the one-source constraints remains unchanged.

Detector | Strategy | livetime | (tracks) | |
---|---|---|---|---|

IceCube | HESE | 1 | 4 yr | |

upgoing | 1 | 6 yr | ||

IceCube-Gen2 | HESE | 10 | 8 yr | |

upgoing | 10 | 12 yr | ||

KM3NeT | HESE | 4 | 8 yr | |

upgoing | 4 | 12 yr |

In future HESE-like analyses, the limits on from the angular power spectrum from IceCube-Gen2 and KM3NeT (summarized in Fig. 5) will outperform point-source searches only if the isotropic flux is due to individually bright sources rarer than . In this hypothetical nondetection scenario, the parameter space associated to blazars would not be constrained much better than it is today using upgoing events (cf. Fig. 4), due to limited improvements in exposure, as well as in angular resolution.

Constraint prospects for future analyses of upgoing (uncontained) tracks are summarized in Fig. 6. In the pessimistic case studied here of a continued nondetection of anisotropy or point sources, KM3NeT and IceCube-Gen2 would (independently and with high significance) rule out a blazar contribution to the high-energy neutrino flux observed today. The angular power spectrum from the next generation of neutrino telescope also has the potential to constrain radio galaxies. Indeed, the upper limits for would reach down to by the time these experiments are decommisioned, well below their neutrino flux anticipated from interactions Hooper2016b . In both HESE and upgoing track analyses, the one-source constraint will still be the most stringent on the population of starburst galaxies, suggesting that it will still be unlikely for the neutrino telescopes to detect them (see also Refs. Feyereisen2016 ; MuraseWaxman ).

These forecast clearly shows that in the future, if the high-energy neutrino sky remains consistent with isotropy, the angular power spectrum will provide much stronger upper limits on the flux of the brightest neutrino source than point-source searches. It also suggests (by comparison with Fig. 4) that if sources are not discovered individually in the near future, they will likely be discovered statistically through the angular power spectrum first. Indeed, due to the respective and scalings of the point-source search and the APS, a statistical discovery becomes increasingly likely the longer point sources are not discovered.

## Vii Conclusions

To conclude, we discussed two constraints on the flux of the brightest neutrino source in the sky, , and how they relate to (or improve on) the null results of the current anisotropy and point-source searches. The one-source limit on manages to reach quite low values, more than one order of magnitude below the existing upper limits based on the search for individual point-like sources in the case of abundant source population such as starburst galaxies. The other approach is based on constraining the variance of source flux distribution (or equivalently, the Poisson angular power spectrum). These upper limits are more powerful for rare source classes, providing complementary information in the case that no source is detected. In particular, analysis of upgoing track events with the current IceCube exposure already has a potential to rule out the scenario of blazar-domination for the diffuse neutrino flux. In addition, the limits based on the angular power spectrum will become more powerful for the next generation of neutrino telescopes. The combination of the two strategies proposed here provide a very efficient way of answering the question: “How bright can the brightest neutrino source be?”

###### Acknowledgements.

We thank Markus Ahlers, John Beacom, Kohta Murase, and an anonymous referee for helpful comments and discussions on the manuscript. This work was supported by the Netherlands Organization for Scientific Research (NWO) through a Vidi Grant.## Appendix A Flux upper limits of the brightest source from point-source searches

The point-source flux upper limits are dependent on declination IceCubePS . In this paper, however, we are interested in a single value of the flux of the brightest neutrino source. Here we shall discuss how we estimate this flux.

Suppose is the flux of the single brightest source somewhere in the sky. Above the flux corresponding to the point-source upper limit at the declination (where ), there will be on average sources from the full sky. The number of sources above this threshold in a declination bin is therefore . We then assign a probability of finding no source brighter than the current point-source upper limits anywhere in the sky, through the Poisson statistics, as

(26) |

By using post-trial 90% CL upper limits from Ref. IceCubePS , , and , we solve this equation for , and obtain .

## Appendix B Relation to source density and luminosity

We shall characterize a source population by its local number density and the neutrino luminosity . Assuming that they are distributed in a local volume where cosmological effects can be neglected, the number of sources that give fluxes greater than is then multiplied by a volume with a radius :

(27) |

from which one can derive . Taking the luminosity distribution into account, we replace with its average over the luminosity function .

Then, as above, the one-source limit is obtained with , which reads

(28) | |||||

Here we again choose .

The break of the flux distribution at its characteristic flux happens when the cosmological expansion comes into play. Although this is dependent on how the source density evolves as a function of redshift and one needs to fully compute in order to be more precise (e.g., Feyereisen2016 ), here we simply approximate that the transition happens at : , where is the luminosity distance. We then obtain using Eq. (7) by replacing measured with , where is a fractional contribution to the measured intensity from the source population. They are

(29) | |||||

(30) |

If the sources are mono-luminous (i.e., the luminosity function is sharply peaked at some value) as is often assumed in the literature Silvestri:2009xb ; Ahlers:2014ioa ; MuraseWaxman , then all these quantities are determined once and are both given. In this case, by equating with the upper limits from the point-source searches, one can place an exclusion line on the plane. In general, however, the luminosity function can range widely, and if it is flatter than , then and hence are mainly sensitive to the upper cutoff of the luminosity function. Such a behavior in the tail region of the luminosity function is typically found for the blazars in the gamma rays BLLacLF ; FSRQLF , and expected in neutrinos too (see the next section).

## Appendix C Examples of blazar models with flat luminosity distribution

If there is a linear correlation between the neutrino and gamma-ray luminosities, , then one can adopt well-established flux distribution from the gamma-ray measurements such as Ref. 2FHL . However, if the neutrinos are produced by the interaction, and if its opacity is dependent of the gamma-ray luminosity, then the scaling can be different from linear. Here, we take recent examples that predict a stronger correlation, , where . There are models of BL Lacs with Tavecchio2015 and flat-spectrum radio quasars (FSRQs) with MID . These kinds of dependence yield a flat distribution of the neutrino luminosities.

The purpose of this section is to obtain the flux distribution starting from the gamma-ray luminosity function, . The neutrino intensity is

(31) |

where is the comoving volume, , and is the luminosity distance corresponding to the redshift . We adopt the luminosity functions from Ref. FSRQLF for FSRQs and Ref. BLLacLF for BL Lacs, but with the cutoff of for the latter case Tavecchio2015 . Using the measured intensity TeV cm s sr IceCubeICRC , we solve Eq. (31) to obtain the constant of proportionality of the scaling relation . Then, the flux distribution is calculated as

(32) |

where both and are now functions of and .

Figure 7 shows for both the models of BL Lac Tavecchio2015 and FSRQs MID , and compare the one of 2FHL 2FHL assuming a linear scaling . All these models are normalized such that each of them can explain the measured diffuse neutrino intensity entirely. This shows that our phenomenological model based on a simple assumption of the broken power law, with at high-flux tail, indeed captures the overall behavior of the flux distribution, predicted with a realistic gamma-ray luminosity function and even in combination with very strong scaling relations between the neutrino and gamma-ray luminosities.

## Appendix D Case of a flat distribution

Here we address the case where and . As seen in the previous section, this case is very difficult to realize, but in order to make our discussion fully generic, we study it. One example of models that can potentially feature a flat tail in the flux distribution is the case where one expects virtually no source in the local volume with . This is again extremely hypothetical and even unrealistic, because even for starburst galaxies, while the redshift evolution is very steep (the luminosity density evolves as or steeper TAM ), the flux distribution has the Euclidean tail, Feyereisen2016 .

In such a case of flat luminosity function exclusively at cosmological distances (), we therefore need to re-derive the relevant equations (7) and (8), as they are based on the assumption of . A schematic representation of the main contributions to the distribution’s first moments is shown in Fig. 8. We find that this time, the contribution to both and is dominated by sources around , and hence, by studying them, we can constrain the flux of the brightest source together. On the other hand, would be entirely unconstrained, even if such a break existed. Also, since the mean intensity is dominated by sources, we do not have to discuss the one-source limit, . Corresponding to Eqs. (7) and (8), we have, for ,

(33) | |||||

(34) |

where , and . Constraints on and are then obtained by solving these equations, given measured and upper limit :

(35) | |||||

(36) |

Again, if a fraction of the total intensity measured is attributed to this source population, then should be replaced with in the equations above.

Figure 9 shows the constraints on and as a function of exposure normalized to that of the 4 years of IceCube operation, for and for HESE events. Rescaling to other values of is trivial by looking at Eq. (35); for and 1.8, we obtain 0.7 and 2 times larger limits on , respectively. Figure 10 is the same as Fig. 9 but for the high-energy upgoing tracks considered in Sec. IV.3, where the exposure is normalized to the current IceCube value with the livetime of 2060 days.

The anisotropy constraints in the case of show that the IceCube neutrinos have to be made by at least tens to hundreds of sources around . The current upper limit on from the angular power spectrum already exceeds the point-source limit. We note that this approach is closely related to a stacking analysis assuming that multiple sources have the same flux, as performed in Ref. IceCubePS . Since the power spectrum is the variance, its sensitivity and hence that to improves linearly with the exposure, while that from the point-source searches goes only as square root of the exposure. This makes the angular power spectrum even more important for the next generation of neutrino telescopes.

## References

- (1) IceCube Collaboration, M. G. Aartsen et al., “Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,” Science 342 (2013) 1242856, arXiv:1311.5238 [astro-ph.HE].
- (2) IceCube Collaboration, M. G. Aartsen et al., “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett. 113 (2014) 101101, arXiv:1405.5303 [astro-ph.HE].
- (3) IceCube Collaboration, M. G. Aartsen et al., “A combined maximum-likelihood analysis of the high-energy astrophysical neutrino flux measured with IceCube,” Astrophys. J. 809 no. 1, (2015) 98, arXiv:1507.03991 [astro-ph.HE].
- (4) IceCube Collaboration, M. G. Aartsen et al., “All-sky search for time-integrated neutrino emission from astrophysical sources with 7 years of IceCube data,” arXiv:1609.04981 [astro-ph.HE].
- (5) IceCube Collaboration, M. G. Aartsen et al., “Searches for small-scale anisotropies from neutrino point sources with three years of IceCube data,” Astropart. Phys. 66 (2015) 39–52, arXiv:1408.0634 [astro-ph.HE].
- (6) IceCube Collaboration, M. G. Aartsen et al., “Observation and Characterization of a Cosmic Muon Neutrino Flux from the Northern Hemisphere using six years of IceCube data,” Astrophys. J. 833 no. 1, (2016) 3, arXiv:1607.08006 [astro-ph.HE].
- (7) IceCube Collaboration, M. G. Aartsen et al., “IceCube-Gen2: A Vision for the Future of Neutrino Astronomy in Antarctica,” arXiv:1412.5106 [astro-ph.HE].
- (8) U. F. Katz, “KM3NeT: Towards a km Mediterranean Neutrino Telescope,” Nucl. Instrum. Meth. A567 (2006) 457–461, arXiv:astro-ph/0606068 [astro-ph].
- (9) M. Ajello et al., “The Cosmic Evolution of Fermi BL Lacertae Objects,” Astrophys. J. 780 (2014) 73, arXiv:1310.0006 [astro-ph.CO].
- (10) M. Ajello et al., “The Luminosity Function of Fermi-detected Flat-Spectrum Radio Quasars,” Astrophys. J. 751 (2012) 108, arXiv:1110.3787 [astro-ph.CO].
- (11) M. Ajello et al., “The Origin of the Extragalactic Gamma-Ray Background and Implications for Dark-Matter Annihilation,” Astrophys. J. 800 no. 2, (2015) L27, arXiv:1501.05301 [astro-ph.HE].
- (12) Y. Inoue and Y. T. Tanaka, “Lower Bound on the Cosmic TeV Gamma-ray Background Radiation,” Astrophys. J. 818 (2016) 187, arXiv:1512.00855 [astro-ph.HE].
- (13) I. Tamborra, S. Ando, and K. Murase, “Star-forming galaxies as the origin of diffuse high-energy backgrounds: Gamma-ray and neutrino connections, and implications for starburst history,” JCAP 1409 (2014) 043, arXiv:1404.1189 [astro-ph.HE].
- (14) M. R. Feyereisen, I. Tamborra, and S. Ando, “One-point fluctuation analysis of the high-energy neutrino sky,” arXiv:1610.01607 [astro-ph.HE].
- (15) M. Di Mauro, F. Calore, F. Donato, M. Ajello, and L. Latronico, “Diffuse -ray emission from misaligned active galactic nuclei,” Astrophys. J. 780 (2014) 161, arXiv:1304.0908 [astro-ph.HE].
- (16) D. Hooper, T. Linden, and A. Lopez, “Radio Galaxies Dominate the High-Energy Diffuse Gamma-Ray Background,” JCAP 1608 no. 08, (2016) 019, arXiv:1604.08505 [astro-ph.HE].
- (17) J. A. Peacock, Cosmological Physics. 1999.
- (18) IceCube Collaboration, M. G. Aartsen et al., “The IceCube Neutrino Observatory - Contributions to ICRC 2015 Part II: Atmospheric and Astrophysical Diffuse Neutrino Searches of All Flavors,” in Proceedings, 34th International Cosmic Ray Conference (ICRC 2015): The Hague, The Netherlands, July 30-August 6, 2015. 2015. arXiv:1510.05223 [astro-ph.HE]. http://inspirehep.net/record/1398539/files/arXiv:1510.05223.pdf.
- (19) M. Ahlers and F. Halzen, “Pinpointing Extragalactic Neutrino Sources in Light of Recent IceCube Observations,” Phys. Rev. D90 no. 4, (2014) 043005, arXiv:1406.2160 [astro-ph.HE].
- (20) K. Bechtol, M. Ahlers, M. Di Mauro, M. Ajello, and J. Vandenbroucke, “Evidence against star-forming galaxies as the dominant source of IceCube neutrinos,” arXiv:1511.00688 [astro-ph.HE].
- (21) IceCube Collaboration, T. Glüsenkamp, “Analysis of the cumulative neutrino flux from Fermi-LAT blazar populations using 3 years of IceCube data,” EPJ Web Conf. 121 (2016) 05006, arXiv:1502.03104 [astro-ph.HE].
- (22) K. Murase, D. Guetta, and M. Ahlers, “Hidden Cosmic-Ray Accelerators as an Origin of TeV-PeV Cosmic Neutrinos,” Phys. Rev. Lett. 116 no. 7, (2016) 071101, arXiv:1509.00805 [astro-ph.HE].
- (23) S. Ando and E. Komatsu, “Anisotropy of the cosmic gamma-ray background from dark matter annihilation,” Phys. Rev. D73 (2006) 023521, arXiv:astro-ph/0512217 [astro-ph].
- (24) S. Ando, E. Komatsu, T. Narumoto, and T. Totani, “Dark matter annihilation or unresolved astrophysical sources? Anisotropy probe of the origin of cosmic gamma-ray background,” Phys. Rev. D75 (2007) 063519, arXiv:astro-ph/0612467 [astro-ph].
- (25) M. Fornasa et al., “The angular power spectrum of the diffuse gamma-ray emission as measured by the Fermi Large Area Telescope and constraints on its Dark Matter interpretation,” Phys. Rev. D94 no. 12, (2016) 123005, arXiv:1608.07289 [astro-ph.HE].
- (26) M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa, “Atmospheric neutrino flux calculation using the NRLMSISE-00 atmospheric model,” Phys. Rev. D92 no. 2, (2015) 023004, arXiv:1502.03916 [astro-ph.HE].
- (27) S. Ando, M. Fornasa, N. Fornengo, M. Regis, and H.-S. Zechlin, “Astrophysical interpretation of the anisotropies in the unresolved gamma-ray background,” arXiv:1701.06988 [astro-ph.HE].
- (28) G. Giacinti, M. Kachelrieß, O. Kalashev, A. Neronov, and D. V. Semikoz, “Unified model for cosmic rays above 10 eV and the diffuse gamma-ray and neutrino backgrounds,” Phys. Rev. D92 no. 8, (2015) 083016, arXiv:1507.07534 [astro-ph.HE].
- (29) C. Righi, F. Tavecchio, and D. Guetta, “High-energy emitting BL Lacs and high-energy neutrinos - Prospects for the direct association with IceCube and KM3NeT,” Astron. Astrophys. 598 (2017) A36, arXiv:1607.08061 [astro-ph.HE].
- (30) P. Padovani, M. Petropoulou, P. Giommi, and E. Resconi, “A simplified view of blazars: the neutrino background,” Mon. Not. Roy. Astron. Soc. 452 no. 2, (2015) 1877–1887, arXiv:1506.09135 [astro-ph.HE].
- (31) P. Padovani, E. Resconi, P. Giommi, B. Arsioli, and Y. L. Chang, “Extreme blazars as counterparts of IceCube astrophysical neutrinos,” Mon. Not. Roy. Astron. Soc. 457 no. 4, (2016) 3582–3592, arXiv:1601.06550 [astro-ph.HE].
- (32) K. Murase, M. Ahlers, and B. C. Lacki, “Testing the Hadronuclear Origin of PeV Neutrinos Observed with IceCube,” Phys. Rev. D88 no. 12, (2013) 121301, arXiv:1306.3417 [astro-ph.HE].
- (33) L. A. Anchordoqui, T. C. Paul, L. H. M. da Silva, D. F. Torres, and B. J. Vlcek, “What IceCube data tell us about neutrino emission from star-forming galaxies (so far),” Phys. Rev. D89 no. 12, (2014) 127304, arXiv:1405.7648 [astro-ph.HE].
- (34) X.-C. Chang and X.-Y. Wang, “The diffuse gamma-ray flux associated with sub-PeV/PeV neutrinos from starburst galaxies,” Astrophys. J. 793 no. 2, (2014) 131, arXiv:1406.1099 [astro-ph.HE].
- (35) N. Senno, P. Mészáros, K. Murase, P. Baerwald, and M. J. Rees, “Extragalactic star-forming galaxies with hypernovae and supernovae as high-energy neutrino and gamma-ray sources: the case of the 10 TeV neutrino data,” Astrophys. J. 806 no. 1, (2015) 24, arXiv:1501.04934 [astro-ph.HE].
- (36) S. Ando, I. Tamborra, and F. Zandanel, “Tomographic Constraints on High-Energy Neutrinos of Hadronuclear Origin,” Phys. Rev. Lett. 115 no. 22, (2015) 221101, arXiv:1509.02444 [astro-ph.HE].
- (37) Fermi-LAT Collaboration, M. Ackermann et al., “Resolving the Extragalactic -Ray Background above 50 GeV with the Fermi Large Area Telescope,” Phys. Rev. Lett. 116 no. 15, (2016) 151105, arXiv:1511.00693 [astro-ph.CO].
- (38) F. Tavecchio and G. Ghisellini, “High-energy cosmic neutrinos from spine-sheath BL Lac jets,” Mon. Not. Roy. Astron. Soc. 451 no. 2, (2015) 1502–1510, arXiv:1411.2783 [astro-ph.HE].
- (39) K. Murase, Y. Inoue, and C. D. Dermer, “Diffuse Neutrino Intensity from the Inner Jets of Active Galactic Nuclei: Impacts of External Photon Fields and the Blazar Sequence,” Phys. Rev. D90 no. 2, (2014) 023007, arXiv:1403.4089 [astro-ph.HE].
- (40) A. Silvestri and S. W. Barwick, “Constraints on Extragalactic Point Source Flux from Diffuse Neutrino Limits,” Phys. Rev. D81 (2010) 023001, arXiv:0908.4266 [astro-ph.HE].
- (41) K. Murase and E. Waxman, “Constraining High-Energy Cosmic Neutrino Sources: Implications and Prospects,” Phys. Rev. D94 no. 10, (2016) 103006, arXiv:1607.01601 [astro-ph.HE].
- (42) Fermi-LAT Collaboration, A. A. Abdo, “Detection of Gamma-Ray Emission from the Starburst Galaxies M82 and NGC 253 with the Large Area Telescope on Fermi,” Astrophys. J. 709 (2010) L152–L157, arXiv:0911.5327 [astro-ph.HE].
- (43) D. Hooper, “A Case for Radio Galaxies as the Sources of IceCube’s Astrophysical Neutrino Flux,” JCAP 1609 no. 09, (2016) 002, arXiv:1605.06504 [astro-ph.HE].