Complementarity of Stacking and Multiplet Constraints on the Blazar Contribution to the Cumulative High-Energy Neutrino Intensity

TL;DR

This study assesses how well blazars can explain the observed high-energy neutrino background by combining gamma-ray luminosity functions, stacking limits, and multiplet constraints, revealing blazars are unlikely to be the dominant sources.

Contribution

It introduces a comprehensive analysis combining stacking and multiplet constraints to evaluate blazar contributions to high-energy neutrinos across various luminosity correlations.

Findings

Blazars are disfavored as the main sources of 100-TeV neutrinos.

Unresolved blazars significantly contribute to the neutrino flux.

Joint constraints limit the neutrino emission models for blazars.

Abstract

We investigate the blazar contribution to the cumulative neutrino intensity assuming a generic relationship between neutrino and gamma-ray luminosities, $L_{\nu}\propto (L_{\rm ph})^{\gamma_{\rm lw}}$. Using the gamma-ray luminosity functions for blazars including flat spectrum radio quasars (FSRQs) and BL Lac objects, as well as the $Fermi$-LAT detection efficiency, we estimate contributions from {blazars resolved by $Fermi$-LAT as well as the unresolved counterpart.} Combining the existing upper limits from stacking analyses, the cumulative neutrino flux from all blazars ({including $Fermi$-LAT resolved and unresolved ones}) are constrained { in the range $0\lesssim\gamma_{\rm lw}\lesssim2.5$}. We also evaluate the effects of the redshift evolution and the effective local number densities for each class of FSRQs, BL Lacs, and all blazars, by which we place another type of constraints on the blazar contribution using the non-detection of high-energy neutrino multiplets. We demonstrate that these two upper limits are complementary, and that the joint consideration of the stacking and multiplet analyses { not only supports the argument that blazars are disfavored as the dominant sources of the 100-TeV neutrino background, but it extends this argument by including also $Fermi$-LAT-unresolved blazars as well, for a more generic luminosity correlation $L_{\nu}\propto (L_{\rm ph})^{\gamma_{\rm lw}}$}.

Equations

$$\frac{d^3{N}_{\mathrm{bl}}}{d{L}_{\mathrm{ph}}dzd\Gamma }={\varphi }_{\mathrm{bl}}(L_{\mathrm{ph}},z)\frac{d{P}_{\mathrm{bl}}}{d\Gamma }\frac{dV}{dz},$$

$$\begin{array}{rl}{\epsilon }_{\nu }{Q}_{{\epsilon }_{\nu }}^{(\mathrm{bl},\mathrm{R})}(z,{\gamma }_{\mathrm{lw}})=& {\int }_{L_{\mathrm{ph},\mathrm{th}}}^{L_{\mathrm{ph},\max }}{\int }_{{\Gamma }_{\min }}^{{\Gamma }_{\max }}{\mathcal{C}}^{-1}{\varphi }_{\mathrm{bl}}(L_{\mathrm{ph}},z)L_{\nu }(L_{\mathrm{ph}})\\ & \times \frac{d{P}_{\mathrm{bl}}}{d\Gamma }d\Gamma d{L}_{\mathrm{ph}}\end{array}$$

$$F_{100}(L_{\mathrm{ph}},z,\Gamma )$$

$$\begin{array}{rl}{\epsilon }_{\nu }{Q}_{{\epsilon }_{\nu }}^{(\mathrm{bl},\mathrm{R})}(z,{\gamma }_{\mathrm{lw}})=& {\int }_{L_{\mathrm{ph},\min }}^{L_{\mathrm{ph},\max }}{\int }_{{\Gamma }_{\min }}^{{\Gamma }_{\max }}{\mathcal{C}}^{-1}{\varphi }_{\mathrm{bl}}(L_{\mathrm{ph}},z)L_{\nu }(L_{\mathrm{ph}})\\ & \times ϵ(L_{\mathrm{ph}},z,\Gamma )\frac{d{P}_{\mathrm{bl}}}{d\Gamma }d\Gamma d{L}_{\mathrm{ph}},\end{array}$$

$$\begin{array}{rl}{\epsilon }_{\nu }{Q}_{{\epsilon }_{\nu }}^{(\mathrm{bl},\mathrm{all})}(z,{\gamma }_{\mathrm{lw}})=& {\int }_{L_{\mathrm{ph},\min }}^{L_{\mathrm{ph},\max }}{\int }_{{\Gamma }_{\min }}^{{\Gamma }_{\max }}{\mathcal{C}}^{-1}{\varphi }_{\mathrm{bl}}(L_{\mathrm{ph}},z)L_{\nu }(L_{\mathrm{ph}})\\ & \times \frac{d{P}_{\mathrm{bl}}}{d\Gamma }d\Gamma d{L}_{\mathrm{ph}}.\end{array}$$

$$E_{\nu }^2{\Phi }_{E_{\nu }}^{(\mathrm{bl},R/all)}({\gamma }_{\mathrm{lw}})=\frac{c}{4\pi }\int dz\frac{{\epsilon }_{\nu }{Q}_{{\epsilon }_{\nu }}^{(\mathrm{bl},R/all)}(z,{\gamma }_{\mathrm{lw}})}{(1+z)}\left|\frac{dt}{dz}\right|,$$

$$\mathcal{F}({\gamma }_{\mathrm{lw}})=\frac{E_{\nu }^2{\Phi }_{E_{\nu }}^{(\mathrm{bl},\mathrm{R})}({\gamma }_{\mathrm{lw}})}{E_{\nu }^2{\Phi }_{E_{\nu }}^{(\mathrm{bl},\mathrm{all})}({\gamma }_{\mathrm{lw}})}.$$

$$E_{\nu }^2{\Phi }_{E_{\nu }}^{(2LAC/3LAC)}=\frac{E_{\nu }^2{\Phi }_{E_{\nu }}^{(2LAC/3LAC,\mathrm{stacking})}}{\mathcal{F}({\gamma }_{\mathrm{lw}})}.$$

$$\begin{array}{rl}{n}_0^{\mathrm{eff}}& to0.0pt\text{lower}4.0pt\sim \text{hss}<1.9\times 10^{-10}\mathrm{Mp}{\mathrm{c}}^{-3}{\left(\frac{{\epsilon }_{\nu }{L}_{{\epsilon }_{\nu }}^{\mathrm{ave}}}{10^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{-3/2}{\left(\frac{b_m{q}_L}{6.6}\right)}^{-1}\\ & \times {\left(\frac{F_{\lim }}{10^{-9.2}\mathrm{GeV}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}}\right)}^{3/2}\left(\frac{2\pi }{\Delta \Omega }\right),\end{array}$$

$$E_{\nu }^2{\Phi }_{E_{\nu }}^{(\mathrm{m})}$$

$${\xi }_z({\gamma }_{\mathrm{lw}})=\frac{\int dz(1+z{)}^{-1}\left|\frac{dt}{dz}\right|f(z,{\gamma }_{\mathrm{lw}})}{\int dz\left|\frac{dt}{dz}\right|},$$

$$n_0^{\mathrm{eff}}=\frac{1}{L_{\nu }^{\mathrm{eff}}}\int d{L}_{\mathrm{ph}}{L}_{\nu }(L_{\mathrm{ph}})\varphi (L_{\mathrm{ph}},0).$$

Full text

1. Introduction

Since the initial detection of high-energy astrophysical neutrinos by the IceCube Neutrino Observatory (Aartsen et al., 2013a, b), a cumulative flux of astrophysical neutrinos in the energy range from $\sim$ 10 TeV to several PeV has been unveiled and measured to a higher precision (Aartsen et al., 2014a, 2015a, 2016). The isotropic distribution of the cumulative flux as well as the background-only results from recent searches for point-like sources and multi-messenger analyses support an extragalactic origin of these neutrinos (Ahlers & Murase, 2014; Aartsen et al., 2014b, 2015b). Up to now, however, the main origin of the cumulative neutrinos still remains unknown.

The flavor ratio measured at Earth, $(\nu_{e}:\nu_{\mu}:\nu_{\tau})\approx(1:1:1)$, is consistent with the prediction from the long-distance oscillations of neutrinos produced through pion decays (Aartsen et al., 2015c), which provides one common framework for the astrophysical models dedicated to explain the cumulative neutrino flux. Many candidates have been proposed and studied (see Ahlers & Halzen, 2015; Mészáros, 2017, for review). Among these candidates, blazars, which are known as a subclass of AGNs with a relativistic jet pointing nearly towards the Earth (Blandford & Rees, 1978; Urry & Padovani, 1995), have been frequently considered as promising ultra-high-energy cosmic-ray (CR) accelerators and high-energy neutrino emitters (e.g., Mücke & Protheroe, 2001; Padovani et al., 2015; Murase, 2017; Resconi et al., 2017). Recently, the IceCube collaboration announced the spatial and temporal coincidence between a muon track neutrino event IceCube170922A and a blazar TXS 0506+056 (IceCube Collaboration et al., 2018a) at the significance $\sim 3\sigma$. Intuitively, if this association is physical, the intimate link between this IceCube neutrino event and the blazar may favor blazars as the main sources of the cumulative neutrino flux, but this is not the case (Murase et al., 2018).

The maximum likelihood stacking searches for cumulative neutrino flux from the second $Fermi$-LAT AGN catalog (2LAC) as well as the point-source searches using the IceCube muon track events and blazars in $Fermi$-LAT 3LAC have independently shown that $Fermi$-LAT-resolved blazars only contribute a small portion of the IceCube cumulative neutrino flux (Aartsen et al., 2017a; Pinat & Sánchez, 2017; Hooper et al., 2018) and the hadronic models of blazar activity are strongly constrained (Neronov et al., 2017), if the specific correlation $L_{\nu}\propto L_{\rm ph}$ is assumed as a prior. Palladino et al. (2019) evaluated the contribution of unresolved sources, and showed that the blazar contribution to the cumulative neutrino flux is constrained unless one makes an ad hoc assumption that lower-luminosity blazars entrain a larger amount of CRs.

Here we argue that, in addition to the stacking analysis, the absence of clustering in high-energy neutrino events, i.e., neutrino multiplets and auto-correleation, can also provide relevant constraints on various classes of proposed sources as the dominant origin of the cumulative neutrino flux (Murase & Waxman, 2016; Ahlers & Halzen, 2014; Aartsen et al., 2014b; Feyereisen et al., 2017; Glauch et al., 2017; Dekker & Ando, 2019). The constraints are sensitive to the redshift evolution of the sources, which are especially powerful for weakly or non-evolving sources such as BL Lac objects (Murase & Waxman, 2016; Murase et al., 2018). But the limits are weaker for rapidly evolving sources such as FSRQs, which could significantly alleviate the constraints, as remarked by Murase et al. (2018). Neronov & Semikoz (2018) studied the constraints on evolving blazar populations and confirmed that fast evolving sources (e.g., $\xi_{z}=5.0$) may indeed relax the neutrino multiplet limits.

In this work, we consider the “joint” implications of these independent analyses for the global blazar population and extend the constraints to a common case where a generic relationship between neutrino and gamma-ray luminosities, e.g., $L_{\nu}\propto(L_{\rm ph})^{\gamma_{\rm lw}}$, is presumed, which is more general than what has been previously considered in such analyses. Physically, the correlation between $L_{\nu}$ and $L_{\rm ph}$ is determined by the interactions between particles and radiation fields inside the sources. Most of physically reasonable models developed on the basis of photohadronic (e.g., $p\gamma$) interactions predict $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$ with indices of 1.0\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}{\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.0 (e.g., Murase et al., 2014; Dermer et al., 2014; Tavecchio & Ghisellini, 2015; Petropoulou et al., 2015; Padovani et al., 2015; Murase & Waxman, 2016; Righi et al., 2017; Murase et al., 2018; Rodrigues et al., 2018). The index ${\gamma_{\rm lw}}$ characterizes the source models and may deviate from this fiducial range for models with increasing complexity. Motivated by this, we treat ${\gamma_{\rm lw}}$ as a free parameter and attempt to reveal the ${\gamma_{\rm lw}}$-dependence of the upper limits on all-blazar contributions. In addition, a new feature of our analysis is that we also consider the effect of $Fermi$-unresolved blazars. One caveat is that, in this study, we assume all sources are equal and emit steadily with a single power-law spectrum. Prior to the the IceCube-170922A alert, IceCube collaboration has found a neutrino excess from the direction of TXS 0506+056 during a 158-day time window in 2014-2015 (IceCube Collaboration et al., 2018b), which reveals the the transient nature of the neutrino emission. We need to keep in mind that the multiplet limits are stronger for flaring sources (Murase et al., 2018). The stacking limits are also applicable to time-averaged emission of the flaring sources, as long as the scaling between neutrino and gamma-ray luminosities hold (Murase et al., 2018).

In the first part (§2), we calculate the ratio of neutrino fluxes from $Fermi$-LAT-resolved blazars and all blazars (including both resolved and unresolved contributions). Combining this ratio with the existing constraints on $Fermi$-LAT-resolved blazars, we estimate the upper limits for all-blazar contributions. The multiplet constraints are given in the second part (§3) where we also derive the effective number densities $n_{0}^{\rm eff}({\gamma_{\rm lw}})$ and the redshift evolution factor $\xi_{z}({\gamma_{\rm lw}})$ for blazars and the subclasses, FSRQs and BL Lacs. In either case, we use the blazar gamma-ray luminosity functions provided by Ajello et al. (2015, 2012, 2014) to reconstruct the neutrino luminosity density. In §4 we conclude with a discussion.

2. Implications of Stacking Limits

Given the differential density of blazars as a function of rest-frame 100 MeV-100 GeV luminosity $L_{\rm ph}$, redshift $z$ and photon index $\Gamma$ defined by the gamma-ray flux $F\propto\varepsilon_{\rm ph}^{-\Gamma}$,

[TABLE]

where the subscript “bl” represents blazars considered in the calculation, $\phi_{\rm bl}(L_{\rm ph},\Gamma)={d^{2}N_{\rm bl}}/{dL_{\rm ph}dV}$ is the luminosity function and $dP_{\rm bl}/d\Gamma$ is the probability distribution of spectral index $\Gamma$, we can directly write down the (differential) luminosity density of neutrinos from $Fermi$-LAT-resolved blazars at redshift $z$,

[TABLE]

where $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$ is the neutrino luminosity, $L_{\rm ph,max}$ is a fixed upper limit of blazar luminosity and the lower limit $L_{\rm ph,th}(L_{\rm ph},z,\Gamma)$ is determined by the $Fermi$ LAT threshold flux $F_{100,\rm th}$ in the energy range of 100 MeV – 100 GeV. In this equation, $\mathcal{C}$ is the normalization coefficient determined by $\varepsilon_{\rm CR,\max}$ and $\varepsilon_{\rm CR,\min}$, the maximum and minimum energy that CRs in blazars can achieve. Since we aim to estimate the neutrino flux from a general luminosity relationship, $L_{\nu}\propto(L_{\rm ph})^{\gamma_{\rm lw}}$, and the physics may be unknown for a general $\gamma_{\rm lw}$, we do not try to provide the details of the gamma-ray and neutrino radiation processes. In this work, we assume that the maximum CR energy is the same for all blazars, as is the normalization factor once the spectral index $s$ of the IceCube neutrino flux is specified.

Here, we present one method to rewrite the integrals in equation 2 by incorporating the $Fermi$-LAT detection efficiency. For a blazar at redshift $z$ with the luminosity $L_{\rm ph}\propto\int_{\varepsilon_{\min}}^{\varepsilon_{\max}}F(\varepsilon)\varepsilon d\varepsilon$, where $\varepsilon_{\max}=100(1+z)\rm\ GeV$ and $\varepsilon_{\min}=100(1+z)\rm\ MeV$, and the photon index $\Gamma$, the integrated photon flux at earth can be written as

[TABLE]

where $d_{L}$ is the luminosity distance between the blazar and the detector. Then the lower limit of the integral in equation 2 can be obtained by requiring $F_{100}(L_{\rm ph,th},z,\Gamma)=F_{100,\rm th}$. Alternatively, thanks to the $Fermi$-LAT detection efficiency $\epsilon(F_{100})$ provided by Abdo et al. (2010), we can simplify equation 2 by using the equivalent detection efficiency $\epsilon(L_{\rm ph},z,\Gamma)=\epsilon(F_{100})$,

[TABLE]

where the lower limit $L_{\rm ph,min}$ reduces to a constant and represents the minimal luminosity of blazars that are considered in this work. To eliminate the instrumental selection effect produced by the low detection efficiency for dimmer blazars and to take all blazars into account, we replace the $L_{\rm ph,th}$ in equation 2 by $L_{\rm ph,min}$, which yields the neutrino luminosity density from all blazars $\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{({\rm bl,all})}(z,{\gamma_{\rm lw}})$, which can be written explicitly as

[TABLE]

Meanwhile, using the LFs for luminosity-dependent density evolution (LDDE) models and parameters provided by Ajello et al. (2012, 2014), we successfully reproduced the redshift evolution of FSRQ and BL Lac luminosity densities illustrated in the Figure 6 of Ajello et al. (2014). At this stage, during the integration of $L_{\rm ph}$, we set the maximum and minimum luminosities to be $10^{50}\ \rm erg\ s^{-1}$ and $10^{40}\ \rm erg\ s^{-1}$, respectively. We also found that the results are consistent with the uncertainties in Ajello et al. (2014) when the limits of the integration were varied by one or two orders of magnitude. Another thing that we need to keep in mind is that we assume the $Fermi$-LAT-unresolved blazars share the identical LFs with the resolved ones. Ackermann et al. (2015) pointed that the index distributions for different blazar classes both for the detected ones and undetected ones are slightly different: the photon spectra of newly-detected FSRQs are slightly softer than the 2LAC ones ($\Delta\Gamma<0.1$) while in contrast there is no significant spectral difference between the two sets of BL Lacs. For the completeness, we also consider a deviation, e.g., 0.2, of the photon spectral index from the best-fit values provided by Ajello et al. (2015, 2012, 2014). Such a test reveals that the resulting $\mathcal{F}(\gamma_{\rm lw})$ remains almost unchanged under a slight derivation of $\Gamma$.

Assuming the neutrino spectra from all blazars have the similar power-law form, e.g., $\varepsilon_{\nu}^{2}\Phi_{\varepsilon_{\nu}}\propto\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{({\rm bl,R/all})}\propto\varepsilon_{\nu}^{2-s}$, and using the comoving neutrino luminosities $\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{({\rm bl,all})}(z,{\gamma_{\rm lw}})$ and $\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{({\rm bl,R})}(z,{\gamma_{\rm lw}})$, the all-flavor neutrino fluxes from $Fermi$-LAT-resolved and all blazars at earth are expected to be

[TABLE]

where $\varepsilon_{\nu}=(1+z)E_{\nu}$. Hence, we can write down the fraction of $Fermi$-LAT-resolved blazars to the cumulative neutrino flux in a simple way that depends only on ${\gamma_{\rm lw}}$,

[TABLE]

Ajello et al. (2015) presented the best-fit parameters in the blazar luminosity functions $\phi_{\rm bl}$, which enables us to compute $\mathcal{F}({\gamma_{\rm lw}})$. Since the redshift correction to the energies leads to one extra term $(1+z)^{2-s}$ to the integrand in equation 6 and another factor $(1+z)^{-1}$ to the integrated flux in equation 3, we conclude that, as a consequence, low-redshift blazars become more important when $s=2.5$, in comparison with the $s=2$ case. Therefore, considering nearby blazars are easier to be detected, a steeper neutrino spectrum predicts a larger $\mathcal{F}({\gamma_{\rm lw}})$, which is confirmed by the thin lines in Figure 1. Moreover, noting that the selection of the minimum and maximum luminosities, e.g., $L_{\rm ph,\min}$ and $L_{\rm ph,\max}$ of a blazar is arbitrary, we tested the reliability of $\mathcal{F}(\gamma_{\rm lw})$ by varying the integral limits and found that the results are not sensitive to $L_{\rm ph,\max}$ and $\mathcal{F}({\gamma_{\rm lw}})$ does not change dramatically in the range {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0 as $L_{\rm ph,min}$ increases from $10^{41}\ \rm erg\ s^{-1}$ to $10^{43}\ \rm erg\ s^{-1}$, as shown in Figure 1. Intuitively, a lower $L_{\rm ph,\min}$ implies that more low-luminosity blazars in the sample are less likely to be detected. Also, for a weaker luminosity dependance ({\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0), the low-luminosity blazars dominate the luminosity density due to the large population. The combined effect is that $\mathcal{F}({\gamma_{\rm lw}})$ decreases in the range {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0. Remarkably, from Figure 1, we can conclude that the contribution from $Fermi$-LAT-resolved blazars is nearly the same as the neutrino flux from all blazars when ${\gamma_{\rm lw}}$ is larger than 1.0. The reason is that, assuming a stronger luminosity dependance ( on other words, a higher $\gamma_{\rm lw}$), the brighter blazars become increasingly important. These high-luminosity blazars have a higher chance to be detected and in this case the neutrinos luminosity densities from $Fermi$-LAT-resolved blazars and all blazars are comparable.

To compute the upper limit of cumulative neutrino flux from all blazars, we use the existing constraints, $E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm 2LAC,stacking)}$ and $E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm 3LAC,stacking)}$, from blazar stacking analyses and point-source searches (Aartsen et al., 2017a; Hooper et al., 2018), which are based on $Fermi$-LAT 2LAC and 3LAC blazars. Combining these existing limits with the fraction of the neutrino flux from $Fermi$-LAT-resolved blazars, we estimate the upper limits of all-blazar contributions from $Fermi$-LAT 2LAC and 3LAC analysis,

[TABLE]

The stacking results themselves have some model dependence. Here, to obtain conservative limits, we adopt the results based on the equal flux weighting for $E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm 2LAC/3LAC,stacking)}$. In general this gives conservative limits, and the luminosity weighting improves the constraints. We will see that, even in this most conservative case, the combined constraints of stacking and multiplet analysis are stringent.

Figure 2 illustrates the upper limits for the all-blazar neutrino flux from $Fermi$-LAT 2LAC and $Fermi$-LAT 3LAC analysis. We show all-flavor neutrino fluxes for all curves and data points in this figure. In the left panel, we assume $s=2$ for the neutrino spectrum. In this case, the stacking analysis of $Fermi$-LAT-2LAC blazars gives 1.2\times 10^{-8}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm 2LAC,stacking)}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.6\times 10^{-8} (in the unit of ${\ \rm GeV\ cm^{-2}\ s^{-1}\ sr^{-1}}$, hereafter). The corresponding upper limits for all blazars calculated using equation 8 are illustrated as the magenta area. The green area in the left panel shows the constraints derived from $Fermi-$LAT 3LAC analysis which predicts 8.0\times 10^{-9}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm 3LAC,stacking)}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.4\times 10^{-8}. For the illustration purpose, we include the IceCube all-flavor neutrino flux 4.8\times 10^{-8}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}E_{\nu}^{2}\Phi_{E_{\nu}}^{(\rm IC)}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}8.4\times 10^{-8} in Figure 2 (the cyan area). To avoid underestimating the upper limits due to the uncertainties of the existing results, we introduced a 50% uncertainty to the constraints derived from stacking analysis, which broadens the areas in the left panel of figure 2. The right panel shows the energy-dependent upper limits for an $\varepsilon_{\nu}^{-2.5}$ neutrino spectrum. The solid lines are obtained by assuming ${\gamma_{\rm lw}}=1.0$ whereas the dashed lines correspond to the case ${\gamma_{\rm lw}}=2.0$. The upper limits from $Fermi$-LAT 2LAC and 3LAC analysis are illustrated as magenta lines and green lines, respectively. In this figure, we showed also the all-flavor neutrinos flux (red points; Aartsen et al., 2015a, 2016), the 6-year high-energy starting events (cyan points; Aartsen et al., 2017c) and the the best fit to the upcoming muon neutrinos scaled to three-flavor case (yellow area). The previous discussion reveals that $\mathcal{F}({\gamma_{\rm lw}})$ may depend on $L_{\rm ph,min}$ moderately, when ${\gamma_{\rm lw}}$ is smaller than 1.0. We will further demonstrate in §3 that, in the range of {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0, the neutrino multiplet constraints are more stringent than the upper limits derived from the stacking analyses, which manifests its complementarity in constraining the cumulative neutrino flux from all blazars over a wide range of $\gamma_{\rm lw}$.

3. Implications of High-Energy Neutrino Multiplet Limits

Here, we present another type of constraints on the origins of IceCube diffuse neutrinos, using the negative results from the clustering test of neutrino-induced muon track events. These high-energy track events are generally detected by IceCube with the angular resolution $\sim 0.5\rm\ deg$, which enables us to determine the incoming directions and perform clustering analysis on their time and spatial distributions. So far, all the clustering tests based on high-energy muon neutrinos have found no statistically significant evidence of clustering in the arrival distribution of neutrinos (Aartsen et al., 2014a; Aartsen et al., 2015d, 2017d, 2019b).

In this section, we investigate the implications of the non-detection of neutrino multiplet sources, and consider the limits on blazar contributions to the cumulative neutrino background. To achieve this goal, we follow Murase & Waxman (2016) and write down the limits on the effective source densities. The formalism presented by Murase & Waxman (2016) is applicable to blazars with a general luminosity weighting $L_{\nu}\propto(L_{\rm ph})^{\gamma_{\rm lw}}$ since the functions $L_{\nu}(dN_{\rm bl}/dL_{\rm ph})\propto(L_{\rm ph})^{\gamma_{\rm lw}+1}\phi_{\rm bl}$ are sharply peaked around some effective luminosities $L_{\rm ph}^{\rm eff}$, which demonstrates that the effective source densities and the neutrino luminosity densities are well defined and constrained. Below, we define these crucial quantities and derive the neutrino multiplet constraints for our blazar case.

Assuming the number of sources that produce more than $k-1$ multiplet events is $N_{m\geq k}$, the constraint from the non-detection of $m\geq k$ multiplet events can be obtained by requiring $N_{m\geq k}\leq 1$. Murase & Waxman (2016) studied the implications to the neutrino sources using the absence of “high-energy” multiplet neutrino sources, and calculated the upper limit on the local source number density for an $\varepsilon_{\nu}^{-2}$ neutrino spectrum,

[TABLE]

where $\varepsilon_{\nu}L_{\varepsilon_{\nu}}^{\rm ave}$ is the time-averaged neutrino luminosity of the source, $F_{\rm lim}\sim(5-6)\times 10^{-10}{\rm\ GeV\ cm^{-2}\ s^{-1}}$ is the 8-year IceCube point-source sensitivity at the 90% confidence level (Aartsen et al., 2017b), $q_{L}\sim 1-3$ denotes a luminosity-dependent correction factor, $\Delta\Omega$ represents the sky coverage of the detector and the details of $m\geq k$ neutrino multiplet constraints are encoded in the factor $b_{m}$. Murase & Waxman (2016) find $b_{m}\simeq 6.6$ for $m\geq 2$ multiplets and $b_{m}\simeq 1.6$ for triplets or higher multiplets (e.g., $m\geq 3$). Note that the point-source sensitivity enters the above expression but the numerical results are obtained by calculating the number of tracks using the muon effective area (Murase & Waxman, 2016).

The purpose of this work is to explore the implications for blazar models using existing equations from previous work without making new analyses on multiplet sources. We simply use the results of the previous analysis by Murase & Waxman (2016), which gives the upper limit on the effective number density, $n_{0}^{\rm eff}(\varepsilon_{\nu}L_{\varepsilon_{\nu}}^{\rm ave})$. Moreover, another reason that we choose this approach is that these results are also consistent with the latest limits on transient sources (after the number density is converted into the rate density, e.g., Aartsen et al., 2019a) 111The limit on the rate density of neutrino transients accounting for the diffuse flux is \rho_{0}^{\rm eff}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1.7\times 10^{4}\,{\rm Gpc^{-3}}\,{\rm yr^{-1}}\,{(b_{m}q_{L}/6.6)}^{2}{({\Delta\Omega/2\pi)}^{2}{(T_{\rm obs}/8~{}{\rm yr})}^{2}}{(\xi_{z}/0.7)}^{-3}\phi^{-3}_{\rm lim,-1}

${\rm max}[N_{\rm fl},1]$, where $N_{\rm fl}\approx f_{\rm fl}T_{\rm obs}/t_{\rm dur}\approx T_{\rm obs}/\Delta T_{\rm fl}$ is the number of flaring periods and $\Delta T_{\rm fl}$ is the typical flare interval (Murase et al., 2018). For \Delta T_{\rm fl}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}T_{\rm obs}, the density and diffuse limits become similar to those for steady sources. Substituting the time-averaged sensitivity gives conservative results because of $F_{\rm lim}>\phi_{\rm lim}/T_{\rm obs}$. For \Delta T_{\rm fl}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}T_{\rm obs}, we expect $\rho_{0}^{\rm eff}T_{\rm obs}\approx n_{0}^{\rm eff}(T_{\rm obs}/\Delta T_{\rm fl})$. Because of n_{0}^{\rm eff}(T_{\rm obs}/\Delta T_{\rm fl})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}n_{0}^{\rm eff}, the limits for steady sources can be regarded as conservative.

One can write the limit on the cumulative neutrino flux from the sources as a function of $n_{0}^{\rm eff}$ and the redshift evolution factor $\xi_{z}$ (Murase et al., 2018):

[TABLE]

where $t_{H}$ is the Hubble time. In this expression, $\xi_{z}$ represents the redshift weighting of the neutrino luminosity of the sources and can be evaluated through (Waxman & Bahcall, 1998)

[TABLE]

where $f(z,{\gamma_{\rm lw}})$ is the redshift evolution function of the neutrino luminosity density normalized to unity at $z=0$ for the luminosity correlation $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$, e.g., for blazars we have $f^{\rm(bl)}(z,{\gamma_{\rm lw}})=[\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{(\rm bl,all)}(z,{\gamma_{\rm lw}})]/[\varepsilon_{\nu}Q_{\varepsilon_{\nu}}^{(\rm bl,all)}(0,{\gamma_{\rm lw}})]$. Similarly, we can also calculate the $\xi_{z}$ for the blazar subclasses, FSRQs and BL Lacs using the luminosity functions from Ajello et al. (2012, 2014). The black, blue and red areas in Figure 3 illustrate the redshift evolution factor $\xi_{z}({\gamma_{\rm lw}})$ for all blazars, FSRQs and BL Lacs, respectively. When ${\gamma_{\rm lw}}=1$, we find $\xi_{z}\sim 7-8$ for the gamma-ray luminosity density evolution of FSRQs and $\xi_{z}\sim 0.6-0.7$ for that of BL Lacs, which are consistent with the values found by Murase et al. (2014) and Murase & Waxman (2016). The solid and dashed boundaries in Figure 3 correspond to the sample schemes, $(L_{\rm ph,min}=10^{42}{\ \rm erg\ s^{-1}},\ L_{\rm ph,max}=10^{50}\ {\rm erg\ s^{-1}})$ and $(L_{\rm ph,min}=10^{40}{\ \rm erg\ s^{-1}},\ L_{\rm ph,max}=10^{52}\ {\rm erg\ s^{-1}})$, respectively. If ${\gamma_{\rm lw}}$ is lower than 1.0, low-luminosity sources at lower redshift contribute a significant component to the total neutrino luminosity density, therefore, a smaller $L_{\rm ph,\min}$ results in a smaller $\xi_{z}$. On the contrary, a strong luminosity correlation with {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1.5 boosts the contribution from high-redshift bright blazars, which leads to a larger $f(z,{\gamma_{\rm lw}})$ at higher redshift and as a result makes $\xi_{z}$ larger, as $L_{\rm ph,max}$ increases.

Besides the factor $\xi_{z}$, it is also necessary to calculate the effective local number density $n_{0}^{\rm eff}$, which characterizes the the number density of sources that dominate the neutrino luminosity density for each specified source population. In this work, we use the luminosity functions in combination with the luminosity weighting relation $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$ to estimate the effective number densities $n_{0}^{\rm eff}$ for blazars, FSRQs and BL Lacs. Here, we follow the procedure presented by Murase & Waxman (2016). For each class of neutrino sources, we define an effective neutrino luminosity $L_{\nu}^{\rm eff}\propto(L_{\rm ph}^{\rm eff})^{\gamma_{\rm lw}}$ using the corresponding effective gamma-ray luminosity $L_{\rm ph}^{\rm eff}$ obtained by maximizing $(L_{\rm ph})^{{\gamma_{\rm lw}}}(dN/d{\ln L_{\rm ph}})=(L_{\rm ph})^{{\gamma_{\rm lw}}+1}\phi(L_{\rm ph},z=0)$, where $\phi(L_{\rm ph},z=0)$ is the local luminosity function of the sources that we are interested in. Since the function $(L_{\rm ph})^{\gamma_{\rm lw}+1}\phi(L_{\rm ph},z=0)$ has a maximum around its extreme point for each source population, we may regard blazars, FSRQs and BL Lacs as “quasi-standard candle” sources, among which the neutrino productions are dominated by the sources distributed closely around one certain effective luminosity $L_{\rm ph}^{\rm eff}$. In this case, we have justified the applicability of the equation appeared in this section to constrain the neutrino fluxes from blazars and the subclasses. The left panel of Figure 4 shows the effective gamma-ray luminosity densities for all blazars (black solid line), FSRQs (blue dashed line) and BL Lacs (red dash-dotted line). Intuitively, $L_{\rm ph}^{\rm eff}$ of FSRQ should be larger than that of BL Lacs since FSRQs are more luminous than BL Lacs. Moreover, the function $(L_{\rm ph})^{{\gamma_{\rm lw}}+1}\phi(L_{\rm ph},z=0)$ achieves its maximum at higher luminosity as ${\gamma_{\rm lw}}$ increases, which naturally explains the monotonic increase of $L_{\rm ph}^{\rm eff}({\gamma_{\rm lw}})$. Considering that low-luminosity BL Lacs dominate the neutrino luminosity density if the luminosity correlation is weak (e.g., {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1) whereas bright FSRQs become increasingly important as ${\gamma_{\rm lw}}$ increases, the blazar effective luminosity $L_{\rm ph}^{\rm eff}$ converges to the BL Lac case when ${\gamma_{\rm lw}}$ is less than 1.0 and then gradually approaches to the FSRQ curve, as is confirmed in Figure 4. With the effective neutrino/gamma-ray luminosity, we can write down the effective local number density of the sources

[TABLE]

The right panel of Figure 4 shows the effective number densities of all blazars (black solid line), FSRQs (blue dashed line) and BL Lacs (red dash-dotted line). As expected, BL Lacs dominate the number density and the blazar effective number density converges to BL Lac and FSRQ curves respectively when {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0 and {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}2.0. Different from $\mathcal{F}({\gamma_{\rm lw}})$ and $\xi_{z}$, $L_{\rm ph}^{\rm eff}$ and $n_{0}^{\rm eff}$ does not depend sensitively on the value of $L_{\rm ph,min}$ and $L_{\rm ph,max}$ in the range 0\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}{\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.5. To interpret this, we need to keep in mind that the former two quantities are determined by the integrations over $L_{\rm ph}$, while $L_{\rm ph}^{\rm eff}$ depends only on the shape/slope of the function $(L_{\rm ph})^{{\gamma_{\rm lw}}+1}\phi(L_{\rm ph},z=0)$. From the left panel of Figure 4, we find that $L_{\rm ph}^{\rm eff}$ lies roughly in the range $10^{43}-10^{49}\rm\ erg\ s^{-1}$ which is covered by the interval $10^{42}-10^{50}{\rm\ erg\ s^{-1}}$, the fiducial range used in our calculation. Meanwhile, the integrand in equation 12 peaks around $L_{\rm ph}^{\rm eff}$, therefore once the peak is included, the effective number density $n_{0}^{\rm eff}$ will not vary too much as the lower and upper bounds of the integral changes.

The above calculations provide the preliminary work and the ingredients needed for calculating the neutrino multiplet limits. Selecting $b_{m}q_{L}\simeq 6.6$ for $m\geq 2$ multiplets and $F_{\rm lim}\simeq 10^{9.2}\ \rm GeV\ cm^{-2}\ s^{-1}$ for an $\varepsilon_{\nu}^{-2}$ neutrino spectrum, the blue dashed, red dashed-dotted and thick black lines in the left panel of Figure 2 illustrate the neutrino multiplet limits for FSRQs, BL Lacs and all blazars, respectively. The blue, red and black areas shows the corresponding uncertainties due to $L_{\rm ph,min}$ and $L_{\rm ph,max}$, as discussed before. From this figure we find that the all-blazar multiplet constraint converges to the FSRQ case at higher ${\gamma_{\rm lw}}$ and to the BL Lac case if ${\gamma_{\rm lw}}$ is less than 1.0, just as expected. We also considered the upper limits for triplet or higher multiplets ($m\geq 3$) by changing the value of $b_{m}q_{L}$ to $1.6$. In this case, the constraints relax to the thin black line. This consequence can be interpreted as the concession of allowing blazars to produce $m=2$ multiplet events. So far, all calculations on the multiplet constraints were based on the $\varepsilon_{\nu}^{-2}$ neutrino spectrum, and to extend the results to a general spectrum, e.g., $s=2.5$, detailed calculations on $F_{\rm lim}$ and $n_{0}^{\rm eff}$ (equation 12) are needed, and our results are conservative in this point. Therefore, in the right panel of Figure 2, only upper limits inferred from stacking analysis are shown.

4. Discussion

In this paper, we considered how two types of analyses, namely stacking and multiplets, constrain on the contribution of blazars to the cumulative neutrino flux, assuming a generalized luminosity weighting $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$. Using the gamma-ray luminosity functions for blazars, FSRQs and BL Lacs, we estimated the ratio of the neutrino fluxes from $Fermi$-LAT-resolved blazars and from all blazars (including unresolved ones), $\mathcal{F}({\gamma_{\rm lw}})$, and the effective number densities, $n_{0}^{\rm eff}({\gamma_{\rm lw}})$, and the redshift evolution factor, $\xi_{z}$, for different source classes. The joint use of a stacking and multiplet analysis, as well as the use of a generalized luminosity function and inclusion of the effect of unresolved blazars, are new aspects which distinguish this analysis from previous ones. The main results are summarized in Figure 2. From this figure we found that the multiplet constraints are the most important at lower values of ${\gamma_{\rm lw}}$, e.g. {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1.0, whereas all-blazar constraints derived from the existing stacking upper limits are more stringent for a stronger luminosity correlation, e.g., \gamma_{\rm lw}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1.5. The joint consideration of these two kinds of limits supports the extended argument that all blazars, including $Fermi-$unresolved ones, are unlikely to dominate the cumulative neutrino background for a generic correlation between the neutrino and gamma-ray luminosities, $L_{\nu}\propto(L_{\rm ph})^{\gamma_{\rm lw}}$, with the index 0\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}\gamma_{\rm lw}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.5. Canonical blazar models, which are physically motivated and based on the leptonic scenario, predict ${\gamma_{\rm lw}}\sim 1.5-2.0$ (Murase et al., 2014). Our results suggest that the stacking constraints are the most stringent for such physically motivated cases. The multiplet and stacking limits are “complementary”, in the sense that these methods have their own advantages in different regimes, and in combination they provide a stronger and tighter constraint than previously, over a wide range of $\gamma_{\rm lw}$, as pointed out by Murase et al. (2018). We also found that while the multiplet constraints are weaker at larger values of ${\gamma_{\rm lw}}$ they become more stringent again for {\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}1.5 due to the rapid decrease of the effective source density.

In this work, we focus on power-law spectra. The limits are stringent for the neutrino flux in the 0.1 PeV range and become weaker at higher energies. For example, neutrino multiplet limits are weaker if one is interested in the origin of $\sim 1$ PeV neutrinos (Murase & Waxman, 2016; Murase et al., 2018; Palladino et al., 2019). It is possible for blazars to explain the dominant fraction of PeV neutrinos by introducing a lower-energy cutoff of the proton maximum energy (Dermer et al., 2014), although neutrinos at 0.1 PeV and lower energies should come from another population of the sources (e.g., Murase et al., 2019).

One of the uncertainties in this work comes primarily from the selection of the lower and upper limits of the luminosity integral, $L_{\rm ph,min}$ and $L_{\rm ph,max}$. As discussed above, we showed that these uncertainties are well controlled, and the final results are reliable if $L_{\rm ph,min}$ and $L_{\rm ph,max}$ are selected in the fiducial ranges $10^{40}-10^{42}\ \rm erg\ s^{-1}$ and $10^{50}-10^{52}\ \rm erg\ s^{-1}$, respectively. From the joint constraints illustrated in Figure 2, we conclude that blazars are disfavored as a dominant source of the cumulative neutrino flux measured by IceCube for a luminosity weighting $L_{\nu}\propto(L_{\rm ph})^{{\gamma_{\rm lw}}}$ with 0.0\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}{\gamma_{\rm lw}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}2.5. Since different blazar models considered for explaining the cumulative neutrino flux can be commonly characterized by the correlation index ${\gamma_{\rm lw}}$ within this range, our calculations on the upper limits and effective number densities would provide rather general constraints for future studies of blazar neutrinos.

We thank Marco Ajello for useful discussion on the usage of the luminosity function and Nick Rodd for the useful communication. The work of K.M. is supported by the Alfred P. Sloan Foundation and NSF grants No. PHY-1620777 and No. AST-1908689, while that of C.C.Y and P.M. is supported by the Eberly Foundation.