Interpretation of the diffuse astrophysical neutrino flux in terms of the blazar sequence
We study if the diffuse astrophysical neutrino flux can come from blazar jets – a subclass of Active Galactic Nuclei (AGNs) – while at the same time respecting the blazar stacking limit based on source catalogs. We compute the neutrino flux from resolved and unresolved sources using an averaged, empirical relationship between electromagnetic spectrum and luminosity, known as the blazar sequence, for two populations of blazars (BL Lacs and FSRQs). Using a source model with realistic neutrino flux computations, we demonstrate that blazars can indeed power the diffuse neutrino flux and obey the stacking limit at the same time, and we derive the conditions for the baryonic loading (proton versus -ray luminosity) evolving over the blazar sequence. Under the hypothesis that blazars power the diffuse astrophysical neutrino flux, we find that the dominant contribution of the diffuse flux up to PeV energies must come from unresolved BL Lacs with baryonic loadings larger than about . A small contribution of resolved high-luminosity FSRQs with baryonic loadings of order can contribute at the highest energies, which is, however, directly tested by the stacking limit. Our results are roughly consistent with the expectations for AGN energetics and have implications for radiation models for AGNs.
The discovery of a diffuse flux of high-energy neutrinos by IceCube Aartsen et al. (2013) has triggered many questions still open nowadays, one of the most relevant being “what is the source of the high-energy neutrinos detected by IceCube?”. Several candidate classes have been proposed to interpret the IceCube data, such as blazars Protheroe (1997); Essey et al. (2010); Murase et al. (2014); Padovani et al. (2016a); Esmaili et al. (2016a), Gamma-Ray Bursts (GRBs) Paczynski and Xu (1994); Vietri (1998); Meszaros and Waxman (2001); Hummer et al. (2012); Murase and Ioka (2013), starburst galaxies Loeb and Waxman (2006); Stecker (2007); Tamborra et al. (2014), and dark matter decay Esmaili and Serpico (2013); Chianese et al. (2016). Some of them are already strongly constrained by measurements, based on the lack of correlations between high-energy neutrinos and known sources. For example, it is known that the contribution of GRBs to the IceCube neutrino signal is of the order of a few percent Aartsen et al. (2017a), whereas the contribution of starburst galaxies may be larger but is also insufficient to explain the observed neutrinos, due to constraints from the extragalactic -ray background Murase et al. (2013); Bechtol et al. (2017).
In this work we study the contribution of blazars to the astrophysical neutrino flux. Blazars are a sub-class of Active Galactic Nuclei (AGNs), which are supermassive black holes that accrete matter from the host galaxy, launching relativistic jets and emitting non-thermal radiation. In the case of blazars, the relativistic jet points in the direction of Earth. Assuming that cosmic rays (CRs) may be accelerated in the jet, we calculate the neutrino emission from CR interactions with the blazar photon fields. To obtain the photon spectrum we rely on the so-called blazar sequence Ghisellini et al. (2017), an empirical relation between the luminosity of a blazar and the photon spectrum emitted.
Considering the full dataset of high-energy starting events (HESE) in IceCube, the contribution of resolved blazars to the astrophysical neutrino flux cannot be larger than Aartsen et al. (2017b) based on the missing association with sources in -ray catalogues. On the other hand, there have been indications of associations of individual neutrino events with AGNs Resconi et al. (2017); Padovani et al. (2016b); Kadler et al. (2016), which means that it is plausible that a few of IceCube’s events stem from resolved blazars. If these indications are confirmed, they will support the AGN origin of the diffuse flux, possibly powered by unresolved objects.
In this work we show that unresolved blazars can provide the dominant contribution to the high-energy neutrinos detected between 100 TeV and a few PeV, while the contribution of very bright sources to the neutrino flux must be highly suppressed in order to respect the blazar stacking limit Aartsen et al. (2017b). This result allows us to draw conclusions about the possible baryonic loading of blazars evolving over the blazar sequence.
In Section II.1 we evaluate the neutrino spectra from blazars numerically, taking into account the differences between the two blazar sub-classes, namely BL Lacs and Flat-Spectrum Radio Quasars (FSRQs). In Sections II.2 and II.3 we discuss the cosmic evolution of BL Lacs Ajello et al. (2014) and FSRQs Ajello et al. (2012), which is necessary to compute the expected diffuse neutrino flux at Earth. We present our results in Section III, where we compute the diffuse flux of neutrinos under three different hypotheses: i) constant baryonic loading, ii) proportionality between the luminosity of blazars in neutrinos () and in -rays (), and iii) a baryonic loading that scales with the source luminosity. A discussion is presented in Section IV, and a conclusion in Section V.
We now describe the methods used to compute the neutrino flux from blazars. In Section II.1 we present the radiation model for the calculation of neutrino spectra across the blazar sequence, in Section II.2 we address the cosmic evolution of blazars, and in Section II.3 we explain the procedure to calculate the diffuse neutrino flux observed at Earth.
ii.1 Neutrinos from the blazar sequence
The calculation of neutrino production in blazars follows closely that in Rodrigues et al. (2018), which is based on the model introduced in Murase et al. (2014). Some details of the present calculation, including key differences to the previous implementation of the model, are discussed in Appendix A. The basic assumption is that CR protons are accelerated in the relativistic jet of the blazar during a period of flaring activity. We assume the same jet speed for all the sources, corresponding to a Lorentz factor . Cosmic rays are injected isotropically during a period of time corresponding to a typical flare, ,111Primed quantities will refer to the rest frame of the jet plasma, while unprimed quantities may refer to the source or the observer’s frame. in a spherical region of size . This is the region represented in dark red in Fig. 1.
Cosmic rays are assumed to undergo second-order Fermi acceleration, with an acceleration rate given by
where is the atomic number of the nucleus being accelerated, is its energy, and is the magnetic field strength in the jet (assumed to scale with the blazar luminosity, see Appendix A). Moreover, is the acceleration efficiency, which is assumed to be in this work. That low value of acceleration efficiency is motivated by matching the primary energy to the cutoff of the IceCube neutrinos, and implies the loss of a possible connection to ultra-high-energy cosmic rays (as e.g. in Rodrigues et al. (2018), where a higher value of was used). We assume CRs are accelerated to a power-law energy distribution in the jet:
where is the maximum energy achieved by the cosmic rays in the source. This is the energy at which cooling or interaction processes become more efficient than acceleration, or the acceleration becomes limited by the size of the blob.
The interaction of CRs with a target photon field is simulated using the NeuCosmA code Baerwald et al. (2012). The main process responsible for neutrino emission is photo-meson production, where in most cases the neutrinos produced carry around 5% of the energy of the primary. The amount of CRs injected in each source is a free parameter of the model; we quantify this quantity by means of the baryonic loading of the source, defined as
where is the total luminosity of injected CRs and is the -ray luminosity of the source, defined here in the range , roughly corresponding to the Fermi-LAT observation range.
The spectral energy distribution (SED) of each source (used as the target photon field for hadronic interactions) is determined according to the most recent implementation of the blazar sequence Ghisellini et al. (2017), an empirical relation based on the Fermi 3LAC blazar catalog Ackermann et al. (2015) that attributes an average SED to each blazar luminosity , distinguishing between BL Lacs and FSRQs. The SEDs are shown in the left panels of Fig. 2: the upper panels refer to BL Lacs, and the bottom panels to FSRQs. All SEDs with and have been extrapolated by renormalizing the brightest and dimmest SEDs (respectively) provided in Ghisellini et al. (2017). As we can see in the top left panel of Fig. 2, high-frequency peaked BL Lacs (HBLs), i.e. BL Lacs with synchrotron and inverse Compton peaks at higher energies, are generally dimmer sources, while low-frequency peaked BL Lacs (LBLs) are brighter. In FSRQs, on the other hand, there is no strong relationship between the frequency of non-thermal emission bands and the jet luminosity (bottom left panel of Fig. 2). Besides synchrotron and inverse Compton emission, FSRQs typically exhibit bright broad lines from atomic emission of the gas surrounding the accretion disk, as well as thermal bumps from an accretion disk and a dusty torus. If we assume that the size of the broad line region (BLR, green region in Fig. 1) increases with luminosity Murase et al. (2014), then the radiation zone in the jet is contained within reach of the external photon fields for high-luminosity (HL) FSRQs (bottom panel of Fig. 1). These external fields, seen in the SEDs of the brightest FSRQs in Fig. 2, enhance CR interactions both inside the jet and during the propagation of escaping CRs through the BLR and near the dusty torus. This is taken into account by means of two additional radiation regions in the modeling of high-luminosity FSRQs (see Appendix A).
In the right panels of Fig. 2 we show the emitted neutrino spectra obtained with our blazar model. The integral of the neutrino spectrum of each blazar yields the total neutrino luminosity of the source, . This quantity is directly proportional to the baryonic loading of the source, since all the emitted neutrinos originate in CR interactions (Fig. 2 refers to the special case ). In general, the emitted neutrino luminosity of a blazar may be written as
where is the product of the neutrino production efficiency and the baryonic loading, defined in Eq. (3). The neutrino production efficiency quantifies the amount of energy from CRs converted into neutrinos in the source (which is independent of the baryonic loading). The quantity will be used throughout this work to express the ratio between the neutrino and -ray luminosity of a given source.
Using the results of Fig. 2, we can calculate the neutrino production efficiency of each source considered in this work. This is shown in the top panel of Fig. 3 as a function of the -ray luminosity of the source. The neutrino production efficiency increases monotonically with the -ray luminosity of the blazar, since higher luminosities imply higher photon densities in the radiation zone in the jet. Note the abrupt increase in the efficiency of FSRQs with , which is due to interactions with external fields, producing additional neutrinos.
In the bottom panel of Fig. 3 we show the maximum energy achieved by protons accelerated in the jet, which is highest for blazars of . Above this luminosity, photo-hadronic interactions of high-energy CRs dominate over acceleration, and the maximum CR energy starts decreasing with luminosity.
ii.2 Source evolution
The calculation of the cumulative neutrino flux from blazars, reported in this work, is based on the distribution of BL Lacs and FSRQs by Ajello et al. Ajello et al. (2014, 2012). In these two works the distribution of blazars is parametrized in order to reproduce the sources observed by Fermi-LAT and the observed diffuse -ray background Abdo et al. (2010). The parametrization is also useful to evaluate the contribution of unresolved sources, i.e. the sources below the Fermi-LAT sensitivity threshold (but expected from a theoretical point of view) that contribute to the diffuse -ray flux.
In Appendix B we give the details of the parameterization by Ajello et al., which we use in Section II.3 to compute the cumulative neutrino flux from blazars. The distribution describes the number of sources per redshift and luminosity and is characterized by 12 parameters (which are different for BL Lacs and FSRQs), reported in Tab. 1.
The distributions of BL Lacs and FSRQs are reported in Fig. 4, in the left and right panels, respectively. In these density plots we represent the quantity , where is the number of sources, is the redshift and . The white curve separates resolved and unresolved sources. It is interesting to notice that FSRQs are almost entirely resolved, since the red region is almost all above the white curve. On the other hand, we notice that with this parameterization there are many BL Lacs expected at high red shift and low luminosity below the white curve, which cannot be detected but can contribute to the total flux of -rays and neutrinos.
In this work we will focus on astrophysical neutrinos, but the distribution by Ajello et al. was originally created to study -ray emission. Using the distribution, it is possible to evaluate the -ray flux expected from BL Lacs and FSRQs. It is interesting to notice that 50% of the -ray flux from BL Lacs is produced by resolved sources and 50% by unresolved sources (see also Tab. 2 of Esmaili et al. (2016b)). For this reason, it is quite natural to expect a similar behavior for neutrinos, and this aspect will be discussed in depth in the next sections. Moreover, if this behavior were confirmed for neutrinos, the lack of correlations between neutrinos and known -ray sources would become less problematic.
ii.3 The expected flux of astrophysical neutrinos
The expected flux of astrophysical neutrinos produced by blazars can be determined using the neutrino flux from a single source (identified by its luminosity and its redshift) and the distribution of blazars in the universe.
The neutrino flux at Earth produced by a single source with -ray luminosity at redshift is given by the following expression:
where is the baryonic loading of the specific source and the neutrino luminosity spectra are represented in Fig. 2 for a baryonic loading . The relation between the total neutrino luminosity of the blazar and its -ray luminosity is given by Eq. (4).
The cumulative flux of neutrinos at Earth is given by convolving the single-source neutrino flux with the distribution of BL Lacs and FSRQs over and redshift , as follows:
where the integration limits and are given in Appendix B. On the other hand, to calculate separately the contribution of resolved and unresolved blazars, the -ray luminosity must be integrated either in the range (unresolved) or (resolved), where is the threshold luminosity for the detectability with Fermi-LAT (corresponding to a flux of Ajello et al. (2017)). The value of depends on the source’s redshift; a good approximation of this function is , with and .
We test three different hypotheses for the baryonic loading and therefore for the relation between the the -ray luminosity and the neutrino luminosity defined in Eq. (4):
- Scenario 1: Constant baryonic loading.
We assume a baryonic loading for all sources.
- Scenario 2: Constant ratio .
We assume that the ratio between neutrino luminosity and -ray luminosity is constant, i.e.. This assumption is model-independent and has been used in previous works, such as Kadler et al. (2016); Righi et al. (2017); Halzen and Kheirandish (2016) to evaluate the flux of neutrinos from BL Lacs. In the context of our model, it implies that the product (cf.Eq. (4)) is constant, which implies that . We allow for different values for BL Lacs and FSRQs, reflecting the potentially different baryonic loadings.
Scenario 3: Baryonic loading varying with as a power law. We assume that the baryonic loading is a power law of , that is universal for BL Lacs and FSRQs. This is a generalization of the second scenario, but it does not allow for different values for BL Lacs and FSRQs.
The theoretical predictions will be compared with the throughgoing muon flux Aartsen et al. (2016), with the high-energy starting events Aartsen et al. (2015) above 100 TeV and with the blazar stacking limit Aartsen et al. (2017b). As discussed in depth in Palladino and Vissani (2016); Palladino et al. (2016); Palladino and Winter (2018), the IceCube data below 100 TeV can be affected by the presence of Galactic neutrinos and atmospheric background (both conventional and prompt neutrinos). For this reason, we choose the throughgoing muons as reference for the extragalactic neutrino flux.
We now discuss the results obtained in the three scenarios described above.
iii.1 Scenario 1: Constant baryonic loading
The diffuse flux obtained choosing a constant baryonic loading is represented in Fig. 5, where a baryonic loading has been chosen in order to match and not overshoot the present IceCube stacking limit concerning blazars Aartsen et al. (2017b). In Fig. 5 the cumulative flux is represented using a purple solid curve, whereas the contributions from BL Lacs and FSRQs are shown using the shaded yellow and cyan regions, respectively. From the same plot we can also see the flux produced by resolved sources (dotted curves), which in this case is not distinguishable from the total flux produced by BL Lacs and FSRQs. This implies that:
Assuming a constant baryonic loading, the flux of neutrinos is fully powered by resolved sources. Therefore, it is not possible to simultaneously interpret the IceCube observations and obey the stacking limit.
A higher baryonic loading would produce a higher neutrino flux, improving the agreement with the observations but violating the stacking limit. Therefore, assuming that astrophysical neutrinos are produced by blazars, a constant baryonic loading is not a viable assumption to interpret the astrophysical neutrinos. Even if one allowed for two different constant values for BL Lacs and FSRQs, the result would not change because the neutrino flux is determined by the resolved sources; therefore, it is not possible to saturate the IceCube measurement obeying the stacking limit at the same time.
iii.2 Scenario 2: Constant ratio
As a second scenario, we consider a constant value for the ratio . Therefore the neutrino luminosity is simply proportional to the -ray luminosity:
We use two different values of for BL Lacs and FSRQs, since these sources have different characteristics that are bound to be reflected in the neutrino production.
As mentioned above, since , the second scenario implies that the baryonic loading scales as
i.e. the baryonic loading decreases with luminosity, compensating exactly for the increase in neutrino production efficiency . We compute the total and resolved fluxes from BL Lacs and FSRQs and we compare them with the measured throughgoing muon flux Aartsen et al. (2016) and with the present stacking limit Aartsen et al. (2017b), for each value of and . The scan of these two parameters is represented in the left panel of Fig. 6, where we show the , , , and regions. The best fit is given by the following set of parameters:
Considering that the efficiencies of neutrino production are comparable for BL Lacs and FSRQs (for objects having a luminosity smaller than ) the meaning of this result is that the baryonic loading of BL Lacs is about 50 times higher than that of FSRQs.
Under the hypothesis of scenario 2, we can interpret the IceCube data without violating the stacking limit, as illustrated in the right panel of Fig. 6. From this plot we notice that the contribution of FSRQs is suppressed (cyan shaded region) while BL Lacs provide the dominant contribution to the diffuse neutrino flux (yellow shaded region). Moreover, the neutrino flux below PeV is dominated by unresolved BL Lacs, whereas the contribution of resolved BL Lacs (dotted yellow curve) emerges at PeV energies and above. Anyway from the left panel of Fig. 6 we notice that, within 1, it is possible to reduce the contribution of BL Lacs to and to increase to contribution of FSRQs to . Using these values the contributions at low energies is still provided by unresolved BL Lacs but the contribution at high energy is provided by FSRQs and no more by resolved BL Lacs.
Since this scenario is able to explain the IceCube data without violating the blazar stacking bound it is natural to analyze a generalization of it. Therefore we relax the strong assumption and we proceed in a more general way in the next section.
iii.3 Scenario 3: Baryonic loading varying with luminosity as a power law
In this section we generalize the previous procedure and assume that the baryonic loading scales with the luminosity as a power law. We define the following function for the baryonic loading, characterized by the presence of two free parameters and :
In this scenario we use only one function for the baryonic loading, describing both BL Lacs and FSRQs. The intrinsic difference between these two classes is contained in , especially for , when FSRQs produce neutrinos efficiently. Moreover it is important to remark that BL Lacs are generally less luminous than FSRQs (as illustrated in Fig. 4); therefore, even using the same function for the baryonic loading of BL Lacs and FSRQs, there will be a difference between the baryonic loading of low-luminosity objects (mostly BL Lacs) and high-luminosity objects (mostly FSRQs).
We do not put any constraint on the parameters of Eq. (10) and we scan them, comparing the theoretical predictions with the measured throughgoing muon flux and the blazar stacking limit, as done in the previous section. The result is shown in the left panel of Fig. 7 and the best fit parameters are obtained by
In this scenario we are also able to explain the throughgoing muon flux without violating the stacking bound. The results and the main conclusion are similar to the ones obtained in the second scenario, with the low-energy part dominated by unresolved sources and the PeV part dominated by resolved sources. The main difference is that in this case the resolved sources are FSRQs and not BL Lacs, but we have to take into account that this consideration is true only for the result obtained considering the best fit parameters, as explained in Section III.2. On the other hand, scenario 2 can be roughly reproduced by using and in scenario 3, although it is not possible to reproduce it exactly since in that case the baryonic loading did not scale as a power law.
In conclusion, the outcome from the second and third scenarios is comparable: the contribution of high-luminosity blazars must be suppressed in order to not violate the stacking limit, whereas the contribution of low-luminosity blazars has to be dominant between roughly 100 TeV and 1 PeV. This essentially means that:
Low-luminosity BL Lacs, such as HBLs, can efficiently produce neutrinos and must have a baryonic loading higher than . On the contrary, FSRQs should have a lower baryonic loading (about less than 100) and therefore they are not efficient in neutrino production.
In order to test the plausibility of our results we check if the baryonic loading that we obtain in the third scenario is compatible with the Eddington luminosity. The Eddington luminosity is the maximum luminosity that a body (such as a star) can achieve when there is balance between the outward radiation force and the inward gravitational force. For accreting black holes of mass , such as AGNs, it is given by:
where is the gravitational constant, is the proton mass and is the Thomson cross section. We choose a Lorentz factor (see Gao et al. (2017)) and , which are typical black hole masses for blazars Ghisellini and Tavecchio (2008); Ghisellini et al. (2010) (see also Fig. 7 of Yu et al. (2015)). Note that the Eddington luminosity is not a hard limit on the expected luminosity of blazars and may be exceeded, especially during flares.
For a jet dominated by baryons, we can use the Eddington luminosity to estimate the maximal plausible baryonic loading by comparing to the physical jet luminosity in -rays (corrected by the beaming factor) as
The comparison between the baryonic loading obtained in Section III.3 and the Eddington luminosity limit is shown in the left panel of Fig. 8, where we represent the the baryonic loading for the best fit in each scenario, and the Eddington limit is represented by the cyan band. From this plot we notice that the baryonic loading obtained is compatible with the Eddington luminosity for the Lorentz factor and the black hole masses that we have discussed. For the sake of completeness, the Eddington limit may vary with , due to the fact that the black hole mass and the ejected luminosity are slightly correlated, as found in Yu et al. (2015). On the other hand, this behavior is contained in the uncertainties that we have used, i.e. within the cyan region in the left panel of Fig. 8; therefore this comparison is sufficient for the purposes of this study.
In the middle panel of Fig. 8 we show the ratio for the three different scenarios. In the case of constant baryonic loading (scenario 1) the ratio (orange curves) is trivial and it is proportional to the neutrino production efficiency (because the baryonic loading is constant). In the second scenario the ratio (red lines) is constant by definition, with two different values for BL Lacs and FSRQs, as discussed in Section III.2. In the last scenario the ratio evolves from approximately 10% to 0.1% for BL Lacs (green curve). FSRQs exhibit a similar behavior up to , where strongly increases due to the increasing of the efficiency in neutrino production.
The right panel of Fig. 8 shows what are the luminosities that contribute more to the neutrino flux. The function is obtained as follows:
where for both BL Lacs and FSRQs. Let us recall from Eq. (4) that , where . The factor refers to the fact that frequency and energy are redshifted. It is clear from this plot that there are many low-luminosity sources that can contribute to the astrophysical neutrino flux and this contribution is significant if the baryonic loading is not constant (both scenarios 2 and 3), as we have found in Sections III.2 and III.3. In any case, we remark that this result is model-dependent since and depend on the specific model for neutrino production, discussed in Section II.1 and Appendix A, and depends on the theoretical cosmic evolution of the sources (which in this work we adopted from Ajello et al. (2014) for BL Lacs and from Ajello et al. (2012) for FSRQs).
In order to test the stability of our result we repeat the calculation of Section III.3 assuming a stacking limit 10 times lower than the present one, since it is plausible that with increasing exposure this limit will become stronger in the future. Under this hypothesis the best fit values for the baryonic loading are:
We present the results in Fig. 9, showing that our model is still satisfying even assuming a much stronger stacking limit. In this situation the contribution of high-luminosity sources must be more suppressed compared to the third scenario (our benchmark), in order to obey the stacking limit. However, the spectral shape of this result does not agree so well with that suggested by the IceCube HESE and the throughgoing muons at PeV energy, since the suppression of the baryonic loading of high-luminosity objects implies the suppression of the flux above PeV energies. This result (together with that obtained in the third scenario considering the present stacking limit) shows that the stacking limit is a powerful tool to extract informations on the baryonic loading of the sources that are supposed to produced high-energy neutrinos.
The results obtained in the second and third scenarios also have effects on the SED modeling, since they suggest that low-luminosity objects are rich in hadrons, which means the second hump of the SED may be produced by decay. For intermediate-luminosity objects the situation is unclear and both hadronic and hybrid (lepto-hadronic) models may be sufficient to explain the -ray and neutrino emission at the same time. High-luminosity objects (such as FSRQs) have a low baryonic loading and, in principle, they could be in agreement with hybrid or even leptonic models. In fact, our model continues to work even if FSRQs do not produce neutrinos at all, as in the second scenario, where the contribution of FSRQs can be null within .
V Summary and conclusion
In this work we have studied the possibility that the diffuse astrophysical neutrino flux is powered by blazars, which are Active Galactic Nuclei viewed in the direction of the jet. A major obstacle has been the AGN stacking limit, which constrains the blazar contribution from the non-observation of correlations between high-energy neutrinos and known (observed) blazars, while unresolved sources may at the same time power most of the diffuse neutrino flux. Using a phenomenological relationship between spectral energy distribution (SED) of blazars and -ray luminosity, known as blazar sequence, and a realistic neutrino production model based on these SEDs, we have derived the implications for blazars assuming that diffuse flux and stacking limit have to be described at the same time.
We have demonstrated that the choice of a constant baryonic loading over the blazar sequence does not allow for a description of the neutrino data, because high-luminosity objects are very efficient neutrino producers; as a consequence, resolved sources will dominate the neutrino flux, and the stacking limit constrains the blazar contribution to the diffuse flux. In order to avoid that, we have allowed the baryonic loading to change as a function of luminosity. We have analyzed two different possibilities: in the first one the ratio between luminosity in neutrinos and luminosity in -rays is constant (implying that the baryonic loading decreases with luminosity), in the second one the baryonic loading scales as a power law of -ray luminosity. We have found that the baryonic loading of low-luminosity objects (such has HBLs, a class of low-luminosity BL Lacs) has to be higher than , whereas the baryonic loading of high-luminosity sources (mostly FSRQs) has to be lower than 100. It is also possible that FSRQs do not contribute to the neutrino flux (within in the scenario ), which would mean these sources are fully leptonic. Low-luminosity objects can then power the diffuse neutrino flux, while the contribution of high-luminosity objects is limited by the stacking limit.
In order to test the plausibility of our results, we have demonstrated that the baryonic loading obtained roughly satisfies the Eddington limit, and we have tested the stability of our result assuming an even stronger stacking limit projected into the future.
In conclusion, we have demonstrated that the observed astrophysical flux of throughgoing muons and the stacking limit for blazars are almost perfectly described if one accepts large enough baryonic loadings for unresolved (low-luminosity) objects, while high-luminosity BL Lacs and FSRQs are disfavored as the main contributors to the diffuse flux of high-energy neutrinos. As a consequence, we expect signatures of hadronic processes in the SED of low-luminosity objects, such as in X-ray or TeV -ray data Gao et al. (2017). This result also alleviates the tension related to the lack of electromagnetic counterparts to the sources of the IceCube neutrinos.
Acknowledgments. A.P. thanks F. Vissani for fruitful past discussions concerning the connection between blazars and high-energy neutrinos. We thank Shan Gao for useful discussions concerning blazar modeling and the distribution of blazars. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 646623).
- Aartsen et al. (2013) M. Aartsen et al. (IceCube), Science 342, 1242856 (2013), eprint 1311.5238.
- Protheroe (1997) R. J. Protheroe, ASP Conf. Ser. 121, 585 (1997), eprint astro-ph/9607165.
- Essey et al. (2010) W. Essey, O. E. Kalashev, A. Kusenko, and J. F. Beacom, Phys. Rev. Lett. 104, 141102 (2010), eprint 0912.3976.
- Murase et al. (2014) K. Murase, Y. Inoue, and C. D. Dermer, Phys. Rev. D90, 023007 (2014), eprint 1403.4089.
- Padovani et al. (2016a) P. Padovani, E. Resconi, P. Giommi, B. Arsioli, and Y. L. Chang, Mon. Not. Roy. Astron. Soc. 457, 3582 (2016a), eprint 1601.06550.
- Esmaili et al. (2016a) A. Esmaili, A. Palladino, and F. Vissani, EPJ Web Conf. 116, 11002 (2016a).
- Paczynski and Xu (1994) B. Paczynski and G. H. Xu, Astrophys. J. 427, 708 (1994).
- Vietri (1998) M. Vietri, Phys. Rev. Lett. 80, 3690 (1998), eprint astro-ph/9802241.
- Meszaros and Waxman (2001) P. Meszaros and E. Waxman, Phys. Rev. Lett. 87, 171102 (2001), eprint astro-ph/0103275.
- Hummer et al. (2012) S. Hummer, P. Baerwald, and W. Winter, Phys. Rev. Lett. 108, 231101 (2012), eprint 1112.1076.
- Murase and Ioka (2013) K. Murase and K. Ioka, Phys. Rev. Lett. 111, 121102 (2013), eprint 1306.2274.
- Loeb and Waxman (2006) A. Loeb and E. Waxman, JCAP 0605, 003 (2006), eprint astro-ph/0601695.
- Stecker (2007) F. W. Stecker, Astropart. Phys. 26, 398 (2007), eprint astro-ph/0607197.
- Tamborra et al. (2014) I. Tamborra, S. Ando, and K. Murase, JCAP 1409, 043 (2014), eprint 1404.1189.
- Esmaili and Serpico (2013) A. Esmaili and P. D. Serpico, JCAP 1311, 054 (2013), eprint 1308.1105.
- Chianese et al. (2016) M. Chianese, G. Miele, S. Morisi, and E. Vitagliano, Phys. Lett. B757, 251 (2016), eprint 1601.02934.
- Aartsen et al. (2017a) M. G. Aartsen et al. (IceCube), Astrophys. J. 843, 112 (2017a), eprint 1702.06868.
- Murase et al. (2013) K. Murase, M. Ahlers, and B. C. Lacki, Phys.Rev. D88, 121301 (2013), eprint 1306.3417.
- Bechtol et al. (2017) K. Bechtol, M. Ahlers, M. Di Mauro, M. Ajello, and J. Vandenbroucke, Astrophys. J. 836, 47 (2017), eprint 1511.00688.
- Ghisellini et al. (2017) G. Ghisellini, C. Righi, L. Costamante, and F. Tavecchio, Mon. Not. Roy. Astron. Soc. 469, 255 (2017), eprint 1702.02571.
- Aartsen et al. (2017b) M. G. Aartsen et al. (IceCube), Astrophys. J. 835, 45 (2017b), eprint 1611.03874.
- Resconi et al. (2017) E. Resconi, S. Coenders, P. Padovani, P. Giommi, and L. Caccianiga, Mon. Not. Roy. Astron. Soc. 468, 597 (2017), eprint 1611.06022.
- Padovani et al. (2016b) P. Padovani, E. Resconi, P. Giommi, B. Arsioli, and Y. L. Chang, Mon. Not. Roy. Astron. Soc. 457, 3582 (2016b), eprint 1601.06550.
- Kadler et al. (2016) M. Kadler et al., Nature Phys. 12, 807 (2016), eprint 1602.02012.
- Ajello et al. (2014) M. Ajello et al., Astrophys. J. 780, 73 (2014), eprint 1310.0006.
- Ajello et al. (2012) M. Ajello et al., Astrophys. J. 751, 108 (2012), eprint 1110.3787.
- Rodrigues et al. (2018) X. Rodrigues, A. Fedynitch, S. Gao, D. Boncioli, and W. Winter, Astrophys. J. 854, 54 (2018), eprint 1711.02091.
- Baerwald et al. (2012) P. Baerwald, S. Hümmer, and W. Winter, Astropart. Phys. 35, 508 (2012), eprint 1107.5583.
- Ackermann et al. (2015) M. Ackermann et al. (Fermi-LAT), Astrophys. J. 810, 14 (2015), eprint 1501.06054.
- Abdo et al. (2010) A. A. Abdo et al. (Fermi-LAT), Astrophys. J. Suppl. 188, 405 (2010), eprint 1002.2280.
- Esmaili et al. (2016b) A. Esmaili, A. Palladino, and F. Vissani, EPJ Web Conf. 116, 11002 (2016b).
- Ajello et al. (2017) M. Ajello et al. (Fermi-LAT), Astrophys. J. Suppl. 232, 18 (2017), eprint 1702.00664.
- Righi et al. (2017) C. Righi, F. Tavecchio, and D. Guetta, Astron. Astrophys. 598, A36 (2017), eprint 1607.08061.
- Halzen and Kheirandish (2016) F. Halzen and A. Kheirandish, Astrophys. J. 831, 12 (2016), eprint 1605.06119.
- Aartsen et al. (2016) M. G. Aartsen et al. (IceCube), Astrophys. J. 833, 3 (2016), eprint 1607.08006.
- Aartsen et al. (2015) M. G. Aartsen et al. (IceCube), Astrophys. J. 809, 98 (2015), eprint 1507.03991.
- Palladino and Vissani (2016) A. Palladino and F. Vissani, Astrophys. J. 826, 185 (2016), eprint 1601.06678.
- Palladino et al. (2016) A. Palladino, M. Spurio, and F. Vissani, JCAP 1612, 045 (2016), eprint 1610.07015.
- Palladino and Winter (2018) A. Palladino and W. Winter (2018), eprint 1801.07277.
- Gao et al. (2017) S. Gao, M. Pohl, and W. Winter, Astrophys. J. 843, 109 (2017), eprint 1610.05306.
- Ghisellini and Tavecchio (2008) G. Ghisellini and F. Tavecchio, mnras 387, 1669 (2008), eprint 0802.1918.
- Ghisellini et al. (2010) G. Ghisellini, F. Tavecchio, L. Foschini, G. Ghirlanda, L. Maraschi, and A. Celotti, mnras 402, 497 (2010), eprint 0909.0932.
- Yu et al. (2015) X. Yu, X. Zhang, H. Zhang, D. Xiong, B. Li, Y. Cha, Y. Chen, X. Huang, and Y. Wang, Astrophys. Space Sci. 357, 14 (2015), eprint 1504.05635.
- Fossati et al. (1998) G. Fossati, L. Maraschi, A. Celotti, A. Comastri, and G. Ghisellini, Mon. Not. Roy. Astron. Soc. 299, 433 (1998), eprint astro-ph/9804103.
- Tavecchio et al. (2010) F. Tavecchio, G. Ghisellini, G. Ghirlanda, L. Foschini, and L. Maraschi, Mon. Not. Roy. Astron. Soc. 401, 1570 (2010), eprint 0909.0651.
- Dermer et al. (2014) C. D. Dermer, M. Cerruti, B. Lott, C. Boisson, and A. Zech, Astrophys. J. 782, 82 (2014), eprint 1304.6680.
- Elvis et al. (1994) M. Elvis, B. J. Wilkes, J. C. McDowell, R. F. Green, J. Bechtold, S. P. Willner, M. S. Oey, E. Polomski, and R. Cutri, Astrophys. J. Sup. 95, 1 (1994).
- Inoue and Tanaka (2016) Y. Inoue and Y. T. Tanaka, Astrophys. J. 828, 13 (2016), eprint 1603.07623.
Appendix A Source model
We present here some details of the source model underlying the present work. While most features are similar to LABEL:~Rodrigues et al. (2018), some parameter are different, which are emphasized in the following discussion. It is also worth noting that in Rodrigues et al. (2018) the model was applied to a previous implementation of the blazar sequence Fossati et al. (1998), based mainly on radio and X-ray observations. The new implementation of the blazar sequence used in this work Ghisellini et al. (2017), on the other hand, was constructed based on the more recent Fermi 3LAC blazar catalog Ackermann et al. (2015), as discussed in Section II.1. Moreover, in Rodrigues et al. (2018) only one sequence is considered, where all low-luminosity sources are BL Lacs and all high-luminosity sources are FSRQs, while in the present work we consider two sequences, one for BL Lacs and another for FSRQs, spanning the same luminosity range.
The magnetic field strength in the jet is considered in this work to be a soft power-law of the blazar luminosity, . This scaling was employed in Appendix A of Rodrigues et al. (2018) because it can yield values of for BL Lacs in the range 0.1-1 G, and for FSRQs in the range to a few G, which are in agreement with previous results (Tavecchio et al. (2010); Dermer et al. (2014)). The scaling is assumed to be the same for BL Lacs and FSRQs, which yields for the brightest FSRQ discussed in this work (see Fig. 2), and for the dimmest BL Lac. In contrast, in the main text of Rodrigues et al. (2018) a stronger scaling was assumed, , which yields much weaker magnetic field strengths in the jet of low-luminosity sources.
In Fig. 10 we explore the photo-hadronic interactions in the jet of two blazars from the blazar sequence; one a high-luminosity FSRQ, and the other a BL Lac of intermediate luminosity (bottom). We compare some of their features with an FSRQ and a BL Lac of similar luminosities, but with the assumptions considered in Rodrigues et al. (2018). In this discussion, as well as throughout this work, we consider the case of CR escape by Bohm-like diffusion, where the rate of escape is proportional to the Larmor radius of the CRs of a given energy. As shown in Rodrigues et al. (2018), the assumption for the particular CR escape mechanism only affects marginally the neutrino spectra ejected by the blazar.
In the left panels of Fig. 10 we show the photon density spectrum in the jet. Note that in the case of the FSRQ (top panel), the external radiation is relativistically boosted into the jet blob, as discussed above (, cf. Fig. 2). In high-luminosity FSRQs (see end of this section), additional components of the SED are considered in addition to non-thermal radiation (assumed to be produced by electrons in the jet). These external components are an infra-red bump from a dusty torus of temperature , an X-ray bump from the corona of the accretion disk Elvis et al. (1994), and two broad lines emitted by molecular gas on the edge of the BLR Murase et al. (2014) (cf. Fig. 1). The broad line emission, as well as a component of the radiation from the dusty torus and the accretion disk, are assumed to isotropize in the BLR Rodrigues et al. (2018), and that component is relativistically boosted into the jet frame. This assumption has also been considered in previous works dealing with hadronic blazar jet models Murase et al. (2014).
In the middle panels of Fig. 10 we show the rates of processes undergone by protons in the jet of the FSRQ (top) and the BL Lac (bottom). The acceleration rate displayed as a solid red line is the one considered in this work, i.e. with an acceleration efficiency of (cf. Eq. (1)). In contrast, we show in dashed red the case of ultra-efficient proton acceleration, as discussed in the main text of Rodrigues et al. (2018). As discussed above, in Rodrigues et al. (2018) low-luminosity jets were considered to have lower magnetic field strengths, which counteracts the effect of a higher acceleration efficiency in the calculation of the acceleration rate (cf. Eq. (1)). The maximum energy achieved by protons in the source (which is proportional to the maximum energy of the emitted neutrinos) is the energy at which acceleration becomes dominated by cooling or hadronic interactions (including adiabatic cooling, synchrotron losses, electron-positron pair production and photo-meson production). The adiabatic cooling timescale is assumed to scale with the size of the region, whose inverse timescale is represented by the gray line. Note that the maximum proton energy marked in the figure corresponds to the acceleration efficiency considered in this work, while it would be 2-3 orders of magnitude higher for .
As we can see by the solid curves in the right panels of Fig. 10, the FSRQ (top) produces neutrinos with a maximum of in the source frame ( in the jet rest frame), and the BL Lac (bottom) produces neutrinos with a maximum of 1 PeV in the source frame. If an acceleration efficiency of were considered (dashed curves), the FSRQ neutrino spectrum would instead peak at 1 EeV and the BL Lac at 50 PeV. This shows that the low acceleration efficiency considered in this work is essential to obtain a neutrino spectrum whose cutoff is compatible with the maximum energy of the observed astrophysical neutrinos. In addition, note that diffusive shock acceleration models Inoue and Tanaka (2016) indicate that high acceleration efficiencies may be difficult to achieve in relativistic shocks, as is the case of blazar jets. On the other hand, with such low acceleration efficiencies, these sources cannot power the diffuse ultra-high-energy CR spectrum (cf. bottom panel of Fig. 3).
Besides the acceleration efficiency, the blob size is another parameter of the model that affects the peak energy and normalization of the ejected neutrino spectrum. While the value used in this work reflects a source with a one-day variability timescale, , other values could be considered. In low-luminosity blazars, the maximum neutrino energy then scales , because cosmic-ray acceleration is limited only by the blob size (see bottom center panel of Fig. 10); in high-luminosity objects, this scaling is weaker because the maximum energy is limited by photo-hadronic interactions (see upper center panel of Fig. 10). The total neutrino luminosity produced by the source scales mainly with the optical thickness to photo-meson production. This grows linearly with and with ; thus, the normalization of the emitted neutrino spectrum scales as .
Finally, the model geometry (briefly discussed in Section II.1) follows the principles discussed in detail in Rodrigues et al. (2018). The jet blob is assumed to lie at a distance to the black hole of (see Fig. 1), which is the same for all blazars. At the same time, the sizes of the BLR and dusty torus ( and ), are assumed to scale with Murase et al. (2014); Rodrigues et al. (2018), where the accretion disk luminosity. In turn, and also scale together as Murase et al. (2014); Rodrigues et al. (2018). Therefore, the jet blob lies inside the BLR only in FSRQs with luminosity . In these cases, we model the BLR and the dusty torus as two additional radiation zones. The CRs that escape the jet then interact with the photons fields present in either zone before leaving the source. For FSRQs with , i.e. with luminosity , radiation from the dusty torus is accounted for in the calculation of CR interactions in the jet, but a one-zone model is employed, such as for BL Lacs and dimmer FSRQs (this is the case for those SEDs in Fig. 2 in which the thermal IR bump, but no optical bumps or emission lines, are present).
Appendix B The distribution of BL Lacs and FSRQs
which represents the number of sources per comoving volume , is the emitted luminosity between 0.1-100 GeV and is the slope of the -ray flux, assuming power law spectra in that energy interval.
For our purpose it is convenient to write the previous expression as a function of the redshift and . The comoving volume is given by:
The comoving distance is defined as follows:
where is the Hubble distance. The function is given by
where , with , . Based on the previous equations, the relation between the comoving volume and the redshift is
The Jacobian obtained from the transformation is given by:
Combining the two previous expressions with Eq. (16) we obtain:
The parametrization by Ajello et al. is characterized by several parameters, which we report in Tab. 1. Only the first two parameters, namely and , are dimensional, respectively in and erg/s. In this discussion, the quantity is always in units of erg/s.
The distribution in the slope is assumed to be Gaussian,
The distribution in luminosity is given by a double power law:
The distribution in is described by an evolution factor, also a double power law:
where and . Summarizing, we have:
which, using Eq. (22), becomes:
The coefficient is a dimensionless number.
Our purpose is to compute the diffuse high-energy neutrino flux from blazars. Since the neutrino spectra at the source only depends on the luminosity, we can integrate the slope (we cannot integrate over the redshift , since the redshift enters in the computation of the neutrino flux at Earth). Henceforth we will use the distribution of BL Lacs and FSRQs defined as follows:
The intervals of integration for BL Lacs are and , and for FSRQs and . The distributions of BL Lacs and FSRQs are represented in Fig. 4.