The following article is Open access

Interpretation of the Diffuse Astrophysical Neutrino Flux in Terms of the Blazar Sequence

, , , and

Published 2019 January 21 © 2019. The American Astronomical Society.
, , Citation Andrea Palladino et al 2019 ApJ 871 41 DOI 10.3847/1538-4357/aaf507

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/871/1/41

Abstract

We study whether the diffuse astrophysical neutrino flux can come from blazar jets—a subclass of active galactic nuclei—while it, at the same time, respects the blazar stacking limit based on source catalogs and is consistent with the observation from TXS 0506+056. 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 at the highest energies 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 low-luminosity 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 105—while only a very small contribution may come from resolved high-luminosity BL Lacs or FSRQs, which can be directly tested by the stacking limit. We find that the blazar TXS 0506+056 is on the verge of these populations in our baseline scenario, at a relatively high luminosity and redshift; as a consequence, we predict about 0.3 γ-ray-neutrino associations per year from the whole population, dominated by BL Lacs with ${L}_{\gamma }\simeq {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and z ∼ 0.1.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The discovery of a diffuse flux of high-energy neutrinos by IceCube (Aartsen et al. 2013) has triggered many questions, 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 predict neutrino fluxes and possibly interpret the IceCube data, such as blazars (Protheroe 1997; Essey et al. 2010; Murase et al. 2014; Esmaili et al. 2016; Padovani et al. 2016; Murase et al. 2018; Righi et al. 2018a, 2018b), gamma-ray bursts (GRBs; Paczynski & Xu 1994; Waxman & Bahcall 1997; Vietri 1998; Meszaros & Waxman 2001; Hummer et al. 2012; Murase & Ioka 2013), starburst galaxies (Loeb & Waxman 2006; Stecker 2007; Tamborra et al. 2014), the cores of active galactic nuclei (AGNs; Stecker 2013), and dark matter decay (Esmaili & 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). A certain contribution to the neutrino flux may also come from the Galactic plane (Palladino & Vissani 2016; Palladino et al. 2016), but the latest measurements of IceCube and ANTARES suggest that this contribution cannot be greater than ∼10% (Albert et al. 2018).

In this work, we study the contribution of blazars to the astrophysical neutrino flux. Blazars are a subclass of AGNs, which are supermassive black holes that accrete matter from the host galaxy, launching relativistic jets and emitting nonthermal 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. As the target particles for CRs, the distributions of those photons in the blazar family are represented by a series of spectra called the blazar sequence (Ghisellini et al. 2017), based on an empirical relation between the luminosity of a blazar and the photon spectrum emitted.

Considering the high-energy starting events (HESEs) in IceCube, the contribution of resolved blazars to the astrophysical neutrino flux cannot be larger than ∼20%–25% (Aartsen et al. 2017b) based on the missing association with sources in γ-ray catalogs, or even less based on theoretical considerations (Palladino & Vissani 2017). On the other hand, there have been indications of associations of individual neutrino events with AGNs (Kadler et al. 2016; Padovani et al. 2016; Resconi et al. 2017), which means that it is plausible that a few of IceCube's events stem from resolved blazars. More recently, direct evidence has been found of a correlation between IceCube neutrinos and the object TXS 0506+056 (Aartsen et al. 2018a, 2018b), which represents a breakthrough discovery in multimessenger astrophysics. This observation may support the AGN origin of the diffuse neutrino flux, possibly powered by unresolved objects.

In this work we show that low-luminosity blazars can provide the dominant contribution to the high-energy neutrinos with energy between a few hundreds of 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. We also discuss explicitly the role of TXS 0506+056 and what we can learn from the population study for future such observations.

In Section 2.1 we evaluate the neutrino spectra from blazars numerically, taking into account the differences between the two blazar subclasses, namely, BL Lacs and flat-spectrum radio quasars (FSRQs). In Section 2.2 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 3, where we compute the diffuse flux of neutrinos under three different hypotheses: (i) constant baryonic loading (defined as the power ratio of injected nonthermal protons to the emitted gamma-rays; see Equation (3)), (ii) proportionality between the luminosity of blazars in neutrinos (Lν) and that in γ-rays (Lγ), and (iii) a baryonic loading that scales with the source luminosity. A discussion is presented in Section 4, focusing on (i) source characteristics predicted by the model as a function of γ-ray luminosity, (ii) implications of the multiplet constraint, i.e., the nonobservation of two neutrino events from the same source, (iii) implications of our model to potential neutrino sources such as TXS 0506+056, and (iv) PKS 0502+049 (He et al. 2018). We conclude the work with Section 5.

2. Methods

We now describe the methods used to compute the neutrino flux from blazars. In Section 2.1 we present the radiation model for the calculation of neutrino spectra across the blazar sequence, in Section 2.2 we address the cosmic evolution of blazars, in Section 2.3 we explain the procedure to calculate the diffuse neutrino flux observed at Earth, and in Section 2.4 we introduce the special blazar TXS 0506+056.

2.1. Neutrinos from the Blazar Sequence

The calculation of neutrino production in blazars follows the same methods used in Rodrigues et al. (2018). 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 Γ = 10. We assume that the source is in a persistent flaring state, corresponding to a duty cycle of 100% (see Section 2.4 for a discussion on the implications of different duty cycle assumptions). CRs are injected isotropically during a period of time corresponding to a typical flare, t' = 106 s,1 in a spherical region of size ${r}_{\mathrm{blob}}^{{\prime} }={{ct}}_{\mathrm{flare}}^{{\prime} }$. This is the region represented in dark red in Figure 1.

Figure 1.

Figure 1. Schematics of the blazar model (not to scale). Top: BL Lac or low-luminosity FSRQ. In these sources the relativistic jet blob, where CRs are injected, lies outside the broad-line region (BLR), and therefore CRs hardly interact with the photon fields produced by the accretion disk or the gas surrounding the BLR—which are deboosted in the blob frame. Bottom: high-luminosity FSRQ (HL-FSRQ), where the acceleration region lies in the BLR and the external photon spectra are relevant for the interactions in the blob frame. If we assume that the luminosity of the source increases with the size of the accretion disk, then in FSRQs brighter than ${L}_{\gamma }=2.4\times {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ the jet blob lies inside the BLR in this model, and the efficiency of photohadronic interactions is increased by the external photon fields of the FSRQ. Figure adopted from Rodrigues et al. (2018).

Standard image High-resolution image

CRs are assumed to undergo second-order Fermi acceleration, with an acceleration rate given by

Equation (1)

where Z is the atomic number of the nucleus being accelerated, E' is its energy, and B' 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 η = 10−3 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 that there is no connection to ultra-high-energy CRs (as, e.g., in Rodrigues et al. 2018, where a higher value of η was used). We assume that CRs are accelerated to a power-law energy distribution in the jet:

Equation (2)

where ${E}_{\max }^{{\prime} }$ is the maximum energy achieved by the CRs 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 photomeson 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

Equation (3)

where LCR is the total luminosity of injected CRs and Lγ is the γ-ray luminosity of the source, defined here in the range 0.1–100 GeV, 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 Lγ, distinguishing between BL Lacs and FSRQs. The SEDs are shown in the left panels of Figure 2: the top panels refer to BL Lacs, and the bottom panels to FSRQs. All SEDs with ${L}_{\gamma }\gt {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and ${L}_{\gamma }\lt {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ 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 Figure 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. For FSRQs, on the other hand, there is no strong relationship between the frequency of nonthermal emission bands and the jet luminosity (bottom left panel of Figure 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 Figure 1) increases with luminosity, then the radiation zone in the jet is contained within reach of the external photon fields for high-luminosity (HL) FSRQs (bottom panel of Figure 1; see Murase et al. 2014). These external fields, seen in the SEDs of the brightest FSRQs in Figure 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).

Figure 2.

Figure 2. Left: SEDs for BL Lacs (top) and FSRQs (bottom) used in this work, given in the rest frame of the jet. The nonthermal SEDs are those from the blazar sequence (Ghisellini et al. 2017), and the external components of the energy spectra are based on Murase et al. (2014). Right: neutrino luminosity spectra calculated for each blazar luminosity, for a baryonic loading of ξ = 1, also in the rest frame of the jet (see Appendix A for details). The label of each curve indicates the logarithm of the γ-ray luminosity in the black hole frame in erg s−1, given by ${L}_{\gamma }={{\rm{\Gamma }}}^{4}\,{L}_{\gamma }^{{\prime} }$, where Γ = 10 is the bulk Lorentz factor of the jet.

Standard image High-resolution image

In the right panels of Figure 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, Lν. This quantity is directly proportional to the baryonic loading ξ of the source, since the emitted neutrinos originate from CR interactions (Figure 2 refers to the special case ξ = 1). In general, the emitted neutrino luminosity Lν of a blazar may be written as

Equation (4)

where ${ \mathcal K }\equiv {\epsilon }_{\nu }\times \xi $ is the product of the neutrino production efficiency and the baryonic loading, defined in Equation (3). The neutrino production efficiency epsilonν ≡ Lν/LCR quantifies the amount of energy from CRs converted into neutrinos in the source (which is independent of the baryonic loading). The quantity ${ \mathcal K }={L}_{\nu }/{L}_{\gamma }$ will be used throughout this work to express the ratio between the neutrino and γ-ray luminosity of a given source.

Using the results of Figure 2, we can calculate the neutrino production efficiency of each source considered in this work. This is shown in the top panel of Figure 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 Lγ ≳ 3 × 1048 erg s−1, which is due to interactions with external fields, producing additional neutrinos.

Figure 3.

Figure 3. Top: neutrino production efficiency of BL Lacs (yellow) and FSRQs (blue) of the blazar sequence as a function of the γ-ray luminosity of the source. Bottom: maximum energy of accelerated protons in the jet of BL Lacs (yellow) and FSRQs (blue) as a function of luminosity. We assume that the radiation zone in the jet is exposed to external photon fields from molecular and thermal emission in the case of FSRQs with luminosity above 3 × 1048 erg s−1 (HL-FSRQs); this is responsible for the jump in neutrino efficiency (see Appendix A).

Standard image High-resolution image

In the bottom panel of Figure 3 we show the maximum energy achieved by protons accelerated in the jet, which is highest for blazars of ${L}_{\gamma }\sim 3\times {10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Above this luminosity, photohadronic interactions of high-energy CRs dominate over acceleration, and the maximum CR energy starts decreasing with luminosity.

2.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 in Ajello et al. (2014, 2012). In these two works the distribution of blazars is parameterized in order to reproduce the sources observed by Fermi-LAT and the observed diffuse γ-ray background (Abdo et al. 2010). The parameterization 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 2.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 Table 2.2

The distributions of BL Lacs and FSRQs are reported in Figure 4. In these figures we represent a set of 9186 sources, obtained from a Monte Carlo extraction in agreement with the distribution presented in Appendix B. The total number of sources is already contained in the parameterization by Ajello et al. We represent the source evolution as a function of redshift z and luminosity ${\ell }\equiv {\mathrm{log}}_{10}[{L}_{\gamma }\,(\mathrm{erg}\,{{\rm{s}}}^{-1})]$. The threshold separating the resolved and unresolved sources is represented by a single value of ${\phi }_{\gamma }=4\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ for simplicity, in order to reproduce the ∼1500 blazars listed in the 3LAC catalog. In reality this value depends strongly on the gamma-ray spectral index, the exposure time, and the position of the source (Ackermann et al. 2015). It is interesting to notice that about 50% of FSRQs are resolved, whereas more than 85% of the expected BL Lacs are not resolved but are expected to contribute to the total flux of γ-rays and neutrinos.

Figure 4.

Figure 4. Left panel: cosmic evolution of BL Lacs and FSRQs as a function of redshift z and logarithm of the γ-ray luminosity Lγ, following the distributions given in Ajello et al. (2014, 2012). The orange curve roughly separates resolved sources (above) from unresolved sources (below). The yellow and cyan numbers denote the number of resolved and unresolved BL Lacs and FSRQs, respectively. Right: distribution of BL Lacs and FSRQs as a function of luminosity, obtained by integrating over the redshift. We notice that FSRQs are accumulated at high luminosity, whereas BL Lacs are characterized by two different populations. The luminosity ${L}_{\gamma }^{\mathrm{evo}}=3.5\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ is indicated.

Standard image High-resolution image

In the right panel of Figure 4 we show the distribution of BL Lacs and FSRQs as a function of their luminosity, obtained integrating on the redshift z. We notice that most FSRQs have a luminosity between ${L}_{\gamma }={10}^{46}$ and 1050 erg s−1, whereas BL Lacs are characterized by two populations with different characteristics: a low-luminosity population (roughly in the range of ${L}_{\gamma }={10}^{44}\mbox{--}{10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$) and a high-luminosity population above ${L}_{\gamma }\simeq {10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. We note that the choice of the upper limit of the integration ${L}_{\gamma }^{\max }={10}^{52}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in Ajello et al. (2014) well exceeds the most luminous blazar observed, but it is consistent with the fact that the region of ${L}_{\gamma }\gt {10}^{50}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ contains no blazars and essentially no contribution to the γ-ray flux, as is shown in Figure 4.

These two populations of BL Lacs are characterized by different cosmological evolutions: the low-luminosity BL Lacs have a negative evolution, i.e., many of these objects are nearby, while high-luminosity BL Lacs have a positive evolution, i.e., are more abundant at high redshifts. Following the parameterization by Ajello et al., we find that the evolution changes from negative to positive when the luminosity is equal to ${L}_{\gamma }^{\mathrm{evo}}\simeq 3.5\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. A summary of these characteristics is reported in Table 1, together with the relative contributions from different blazar classes to the resolved and unresolved γ-ray flux.

Table 1.  Summary of Some Characteristics of BL Lacs and FSRQs: The Evolution, the Number of Resolved and Unresolved BL Lacs and FSRQs, and the Resolved vs. Unresolved γ-ray Flux Contribution, Listed Separately for Each Population

  Evolution No. Resolved Sources No. Unresolved Sources Resolved Flux Unresolved Flux
Low-luminosity BL Lacs Negative 359 6070 64% 36%
High-luminosity BL Lacs Positive 609 981 90% 10%
FSRQs Positive 566 601 97% 3%
All blazars 1534 7652 88% 12%

Download table as:  ASCIITypeset image

In this work we focus on astrophysical neutrinos, whereas the distribution by Ajello et al. was originally created to study the γ-ray emission. Indeed, using this distribution, it is possible to evaluate the γ-ray flux expected from BL Lacs and FSRQs. As stated before, a non-negligible fraction of this γ-ray flux can be produced by unresolved sources. Therefore, 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 the contribution of unresolved sources to the neutrino flux were large enough, the lack of correlations between neutrinos and known γ-ray sources would become less problematic.

2.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 Φs produced by a single source with γ-ray luminosity Lγ at redshift z is given by the following expression:

Equation (5)

where ξ() is the baryonic loading of the specific source, the neutrino luminosity spectra EdLν/dE are represented in Figure 2 for a baryonic loading ξ = 1, and Dc(z) is the comoving distance, defined in Appendix B. The relation between the total neutrino luminosity of the blazar and its γ-ray luminosity is given by Equation (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 Lγ and redshift z, as follows:

Equation (6)

where the integration limits 1, 2, z1, and z2 are given in Appendix B. On the other hand, to calculate the contribution of resolved and unresolved blazars separately, the γ-ray luminosity must be integrated either in the range [1, vis] (unresolved) or in the range [vis, 2] (resolved), where vis is the average luminosity threshold discussed in Section 2.2, corresponding to a flux of ${\phi }_{\gamma }^{\mathrm{vis}}=4\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$. Therefore, vis is a function of the redshift, namely, ${{\ell }}_{\mathrm{vis}}{(z)=4\pi \times {\phi }_{\gamma }^{\mathrm{vis}}\times {D}_{c}(z)}^{2}{(1+z)}^{2}$.

It is important to remark that the previous procedure implicitly contains the hypothesis that a source emits always an average luminosity Lγ.

2.4. Neutrinos from the Blazar TXS 0506+056

In 2017 September a high-energy neutrino, event IceCube-170922A, was detected by IceCube (Aartsen et al. 2018b). Based only on its energy and direction, this event has a probability of 50% to be an astrophysical neutrino rather than having an atmospheric origin. On the other hand, after the neutrino detection different γ-ray telescopes have detected a flaring state of the source TXS 0506+056 (Aartsen et al. 2018a), which lies in a direction consistent with that of the IceCube event. This correlation in space and time increased the probability of this event being of astrophysical origin, excluding the atmospheric hypothesis at the ∼3.5σ level.

The aforementioned astrophysical source is considered to be a BL Lac, located at redshift zTXS = 0.3365 ± 0.0010, with a luminosity of ${L}_{\gamma }\simeq 1.3\times {10}^{47}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ during the weeks before and after the neutrino event (Aartsen et al. 2018a). The luminosity averaged over all historical data is ${L}_{\gamma }\simeq 2.8\times {10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Aartsen et al. 2018a), as marked in the right panel of Figure 4. From here on we will use this time-averaged luminosity to study the characteristics of TXS 0506+056 within the blazar population.

This discovery represents a breakthrough in the field of multimessenger astronomy, since it is the first evidence of correlation between high-energy neutrinos and γ-rays. Several theoretical papers have been written after this discovery attempting to explain the neutrino emission with blazar radiation models (Ahnen et al. 2018; Cerruti et al. 2018; Gao et al. 2018; Keivani et al. 2018; Liao et al. 2018; Padovani et al. 2018). It is therefore important to analyze the role of this source in the context of our general description, which focuses on the entire blazar population. Note that, looking back at historical data, IceCube has also found an excess from the same direction between 2014 September and 2015 March, with respect to atmospheric backgrounds. This constitutes 3.5σ evidence of neutrino emission from the direction of TXS 0506+056, independent of and prior to the 2017 flaring episode (Aartsen et al. 2018b). On the other hand, the FSRQ PKS 0502+049 displayed flaring activity around the time of the neutrino flare observed by IceCube between the 2014 and 2015 flares, and its position is compatible with the angular resolution of the events (He et al. 2018; Padovani et al. 2018), which makes this source another interesting test case for our model. In Sections 4.3 and 4.4 we discuss the implications of the current work for both these sources.

3. Results

In this section we analyze the consequences of the hypothesis that the most energetic neutrinos observed by IceCube are produced by blazars. This means that the neutrino flux at the highest energies is saturated by blazars, while at lower energies additional contributions may be present owing to the expected peaked nature of the blazar neutrino flux. Although other studies (Aartsen et al. 2017b; Hooper et al. 2018; Murase et al. 2018) disfavor blazars as the dominant source of the astrophysical neutrinos based on different arguments, it is still premature to rule out this source class in light of the potential contribution from unresolved blazars and the fact that they may dominate only at the highest energies. One motivation for such an investigation is that blazars dominate the gamma-ray sky above 100 GeV; if these gamma rays are of hadronic origin, the same sources ought to efficiently produce neutrinos.

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 Equation (4):

  • Scenario 1: Constant baryonic loading: We assume a constant baryonic loading $\xi ({\ell })=\tilde{\xi }$ for all sources, as in Zhang & Li (2017).
  • Scenario 2: Constant ratio ${L}_{\nu }/{L}_{\gamma }$: We assume that the ratio between neutrino luminosity and γ-ray luminosity is constant, i.e., ${L}_{\nu }/{L}_{\gamma }\equiv { \mathcal K }=\mathrm{const}$. This assumption is model independent and has been used in previous works, such as Kadler et al. (2016), Righi et al. (2017), Halzen & Kheirandish (2016), and Wang & Li (2016), to evaluate the flux of neutrinos from BL Lacs. In the context of our model, it implies that the product $\xi \times {\epsilon }_{\nu }={ \mathcal K }$ (see Equation (4)) is constant, which means that $\xi \propto {({\epsilon }_{\nu })}^{-1}$. We allow for different values for low-luminosity BL Lacs and high-luminosity BL Lacs/FSRQs (for the purposes of this analysis we consider high-luminosity sources as one group), reflecting the potentially different baryonic loadings.
  • Scenario 3: Baryonic loading evolving with ${L}_{\gamma }$: We assume that the baryonic loading changes continuously as a function of ; this function is assumed to be universal for BL Lacs and FSRQs. This is a generalization of the second scenario.

The theoretical predictions will be compared with the throughgoing muon flux (Aartsen et al. 2016), with the HESEs (Aartsen et al. 2015) above 100 TeV, and with the blazar stacking limit (Aartsen et al. 2017b). As discussed in depth in Palladino & Vissani (2016), Palladino et al. (2016), and Palladino & Winter (2018), the IceCube data below 100 TeV can be affected by the presence of Galactic neutrinos and residual 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.

3.1. Scenario 1: Constant Baryonic Loading

The diffuse flux obtained choosing a constant baryonic loading is represented in Figure 5, where a baryonic loading $\tilde{\xi }=150$ has been chosen in order to match and not overshoot the present IceCube stacking limit for blazars (Aartsen et al. 2017b). Let us remark that the stacking limit has been obtained by IceCube assuming a power-law spectrum, although we expect a different spectrum for blazars. A new preliminary analysis has been presented in the conference proceedings (Aartsen et al. 2017c) based on different models. However, these models are not consistent with each other and yield very different predictions. We therefore consider in this work the only published stacking limit, choosing that obtained for a spectral index of 2.2, which is roughly that of the IceCube throughgoing muons. In Figure 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 (dashed curves), which in this case is not distinguishable from the total flux produced by BL Lacs and FSRQs. This implies the following:

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.

Figure 5.

Figure 5. Scenario 1: constant baryonic loading for all sources. The black points represent the IceCube HESEs (all flavor; Aartsen et al. 2015), the green line represents the flux of throughgoing muons (x3; Aartsen et al. 2016), and the red line represents the current IceCube blazar stacking limit (Aartsen et al. 2017b). The shaded yellow and cyan regions denote the contribution of BL Lacs and FSRQs, respectively, to the total flux of neutrinos, whereas the dotted yellow and cyan curves denote the contribution from resolved sources only. The purple solid curve represents the total neutrino flux expected from blazars. In this case most of the flux is powered by resolved sources, particularly FSRQs.

Standard image High-resolution image

3.2. Scenario 2: Constant Ratio ${L}_{\nu }/{L}_{\gamma }$

As a second scenario, we consider a constant ratio ${ \mathcal K }\equiv {L}_{\nu }/{L}_{\gamma }$. We assume two different values for low-luminosity BL Lacs (${{ \mathcal K }}_{\mathrm{low} \mbox{-} \mathrm{lum}}$) and high-luminosity BL Lacs/FSRQs (${{ \mathcal K }}_{\mathrm{high} \mbox{-} \mathrm{lum}}$). This assumption comes from the following consideration: in order to power the neutrino flux detected by IceCube without violating the stacking limit, the contribution of unresolved sources (mainly low-luminosity sources) has to be enhanced, and, at the same time, the contribution of resolved sources (mainly high-luminosity sources) has to be suppressed. As a splitting point, we choose ${L}_{\gamma }^{\mathrm{evo}}$, which possibly separates (for BL Lacs) two different populations; see the right panel of Figure 4 and Section 2.2. Therefore, we have

where ${L}_{\gamma }^{\mathrm{evo}}=3.2\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$.

Since in scenario 2 the value of ${ \mathcal K }={\epsilon }_{\nu }\times \xi $ is constant for all sources, for our source model this implies that the baryonic loading scales as $\xi {({\ell })\propto ({\epsilon }_{\nu }({\ell }))}^{-1}$, i.e., the baryonic loading ξ() decreases with luminosity in a way that it exactly compensates for the increasing neutrino production efficiency epsilonν(). Note that the behavior of the baryonic loading depends on the particular model of neutrino production in the source, while the assumption Lν ∝ Lγ is not model dependent.

We compute the total and resolved fluxes from BL Lacs and FSRQs and 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 ${{ \mathcal K }}_{\mathrm{low} \mbox{-} \mathrm{lum}}$ and ${{ \mathcal K }}_{\mathrm{high} \mbox{-} \mathrm{lum}}$. The scan of these two parameters is represented in the left panel of Figure 6, where we show the 1σ, 2σ, 3σ, and ≥3σ regions. The best fit is given by the following set of parameters:

Equation (7)

The interpretation of that result is that the neutrino flux must be powered by low-luminosity sources in order to avoid overshooting the stacking limit. Within 1σ, ${{ \mathcal K }}_{\mathrm{high} \mbox{-} \mathrm{lum}}\leqslant 0.5 \% $ is allowed, providing an upper limit to the contribution of high-luminosity sources.

Figure 6.

Figure 6. Scenario 2: constant ratio ${ \mathcal K }={L}_{\nu }/{L}_{\gamma }={\epsilon }_{\nu }\times \xi $. Left: scan of the values of ${ \mathcal K }={\epsilon }_{\nu }\times \xi $ for low- and high-luminosity sources, taking into account both IceCube observations and the stacking limit. The number of standard deviations σ is calculated for 2 degrees of freedom. The best fit is marked with a black star. Middle: theoretical expectation for the neutrino flux, as described in Figure 5, compared to the IceCube results. The flux of neutrinos is obtained considering the best-fit values of epsilonν × ξ, i.e., epsilonν × ξ = 10.5% for low-luminosity sources and epsilonν × ξ = 0% for high-luminosity sources. In this scenario, the neutrino flux is powered by both unresolved and resolved sources. Right: same as the middle panel, including the effect of other possible contributions (multicomponent model), such as the residual atmospheric background and Galactic neutrinos, as discussed in Palladino & Winter (2018).

Standard image High-resolution image

Within scenario 2, we can interpret the IceCube data at the highest energies without violating the stacking limit, as illustrated in the middle panel of Figure 6. From this plot we notice that the contribution of FSRQs is absent while BL Lacs (particularly those with low luminosity) provide the dominant contribution to the diffuse neutrino flux (yellow shaded region). Combining this result with the idea of the multiple components presented in Palladino & Winter (2018), which may power the neutrino flux at lower energies, we obtain the total neutrino flux (shown in the right panel of Figure 6), as a combination of (i) the neutrino flux from blazars and (ii) the neutrino flux coming from residual background (atmospheric neutrinos, ϕ(E) ∼ E−3.7 below 100 TeV) and Galactic neutrinos (ϕ(E) ∼ E−2.6 with energy cutoff at 150 TeV); we denote these additional neutrinos as "other components." The additional component has not been refitted in this paper, but it is exactly the same as reported in Palladino & Winter (2018). It is important to remark that in principle the atmospheric background should have been already subtracted in the IceCube data points. On the other hand, in Palladino & Winter (2018) we discuss the possibility that a certain residual background can still affect the IceCube measurement, even after the subtraction of the background (see particularly Figure 2 of Palladino & Winter 2018). The residual atmospheric background may become relevant especially below 100 TeV in the HESE analysis. We notice from the right panel of Figure 6 that the combination of scenario 2 and these additional components allows for an interpretation of the shape of the spectrum measured by IceCube.

The result of scenario 2 is much more optimistic than that of Hooper et al. (2018). The main reason is that Hooper et al. (2018) only consider the contribution of detected sources, presenting an updated version of the analysis contained in Aartsen et al. (2017b). On the contrary, we have shown that under the hypotheses of scenario 2, the contribution of unresolved sources can be relevant.

3.3. Scenario 3: Baryonic Loading Evolving with Luminosity as a Power Law

Here we generalize the previous scenario and assume that the baryonic loading scales with the luminosity as a continuous function, defined as follows:

  • 1.  
    For low-luminosity sources, we roughly replicate ${{ \mathcal K }}_{\mathrm{low} \mbox{-} \mathrm{lum}}\equiv {L}_{\nu }/{L}_{\gamma }=10.5 \% $ suggested by the best-fit result of scenario 2. Therefore, the baryonic loading is given by $\xi ({\ell })\simeq {{ \mathcal K }}_{\mathrm{low} \mbox{-} \mathrm{lum}}/{\epsilon }_{\nu }({\ell })$ for our production model.
  • 2.  
    For high-luminosity sources we use the information ${{ \mathcal K }}_{\mathrm{high} \mbox{-} \mathrm{lum}}\equiv {L}_{\nu }/{L}_{\gamma }\lt 0.5 \% $ from scenario 2, and we derive an upper limit for the baryonic loading, indicated above.
  • 3.  
    For the region in between, we use the information from TXS 0506+056 in the model presented in Gao et al. (2018), where the baryonic loading of the source TXS 0506+056 is calculated to be ξ ≃ 3 × 104 during the flare.

The information provided above can be checked looking at the middle panel of Figure 8, where the function ${ \mathcal K }$ is represented as a function of luminosity. Let us note that at low luminosity the value of ${ \mathcal K }$ for FSRQs is not equal to ≃10%, as stated before. On the other hand, we have to take into account that at low luminosity there are no FSRQs at all (see the right panel of Figure 4); therefore, this fact cannot produce any effect on the calculation.

The inclusion of the information from TXS 0506+056 is relevant for our purpose, since TXS 0506+056 has been identified as the first possible source of an IceCube neutrino event (Aartsen et al. 2018a, 2018b). Although it may not be representative for the whole population, it is so far the only direct evidence of a correlation between neutrinos and γ-rays, and therefore it is the best available piece of information at this point. As a consequence, our model is consistent with TXS 0506+056 by construction.

As in scenario 2, we can describe 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, as can be observed in Figure 7. However, this scenario allows for a small contribution by FSRQs, which was not present in the last scenario. The actual function obtained for the baryonic loading is represented by the blue curve in the left panel of Figure 8. On the other hand, in this case we have chosen the upper limit (within 1σ) of the baryonic loading for high-luminosity sources; therefore, the contribution of FSRQs can also be lower or even negligible, as can be seen by the blue shaded region in the left panel of Figure 8, representing the uncertainty in the baryonic loading. Note that TXS 0506+056 lies on the verge between low-luminosity sources, which dominate the diffuse neutrino flux, and the flux-limited high-luminosity sources; it is therefore a rather special case in the context of our model.

Figure 7.

Figure 7. Scenario 3: baryonic loading changing continuously with luminosity, evolving from a baryonic loading larger than 105 for low-luminosity sources to a baryonic loading smaller than 10 for high-luminosity sources; see main text. As in scenario 2, the neutrino flux is powered by both unresolved and resolved sources, with a dominant contribution from low-luminosity BL Lacs.

Standard image High-resolution image
Figure 8.

Figure 8. Left: baryonic loading as a function of Lγ (best fit as blue curve; uncertainties represented by blue band). The region disfavored by the Eddington luminosity is represented in gray. The red region represents the Eddington luminosity according to a black hole mass range ${10}^{7.5}$${10}^{9.5}\,{M}_{\odot }$ with a Lorenz factor range of Γ = 10–40. Middle: neutrino to gamma-ray luminosity as a function of Lγ. Right: same quantity as in the middle panel multiplied by dN/dℓ, the differential number of sources of luminosity . All panels are given for scenario 3.

Standard image High-resolution image

In conclusion, the results from the second and third scenarios are 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. This conclusively means the following:

Low-luminosity BL Lacs can be bright neutrino sources, and they should have a baryonic loading higher than 105 in order to power the neutrino events in IceCube. On the contrary, FSRQs must have a much lower baryonic loading (about 10) and are therefore dim in the neutrino channel.

From here on this scenario will represent our baseline model, and its implications will be discussed in the next section.

4. Discussion

In this section we discuss some aspects to test the plausibility of our baseline model (scenario 3), namely, (i) the evolution with γ-ray luminosity of the baryonic loading and neutrino luminosity of blazars, including a comparison of the baryonic loading with the Eddington luminosity; (ii) the discussion concerning the multiplet constraint, i.e., the absence of two neutrinos coming from the same source; (iii) the comparison between our general picture and the expectations from the blazar TXS 0506+056, which arouses at present great interest owing to the recent evidence for a correlation between one of its flares and one IceCube neutrino detected in 2017 September; and (iv) a discussion concerning the object PKS 0502+049 as a neutrino emitter.

4.1. Source Parameters with γ-ray Luminosity

In order to test the plausibility of our result, we check whether the baryonic loading that we obtain 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, such as AGNs, of mass M, it is given by

Equation (8)

where G is the gravitational constant, mp is the proton mass, and σT is the Thomson cross section. For the Lorentz factor Γ, a large range of values has been used in the literature of blazar modeling (e.g., Boettcher et al. 2013). Here we have chosen a conservative range between 10 and 40. A range of supermassive black hole masses of $M\in [{10}^{7.5},{10}^{9.5}]\,{M}_{\odot }$ is assumed, which is typical for blazars (Ghisellini & Tavecchio 2008; Ghisellini et al. 2010; see also Figure 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 (Sadowski & Narayan 2015).

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

Equation (9)

The comparison between the baryonic loading ξ() obtained in Section 3.3 and the Eddington luminosity limit is shown in the left panel of Figure 8, where both the best fit (blue curve) and the uncertainties (blue band) are represented. 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. Note that the quantity Ledd may also vary with Lγ, 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 red region in the left panel of Figure 8; therefore, this comparison is sufficient for the purposes of this study. It is, however, important to remark that Equation (9) is a necessary, rather than sufficient, condition. For example, in unified schemes of radio-loud AGNs (Urry & Padovani 1995), in which BL Lacs are associated with FR I radio galaxies with radiatively inefficient accretion flows that do not produce broad emission lines, we expect highly sub-Eddington accretion powers and similarly sub-Eddington jet powers (Sikora et al. 2007). While these studies do not seem to be compatible with our findings, one may speculate that the baryonic injection comes during flares only, when the Eddington luminosity may in fact be exceeded (Gao et al. 2018).

In the middle panel of Figure 8 we show the ratio Lν/Lγ for our baseline scenario 3. The ratio Lν/Lγ evolves from approximately 10% to 0.01% for BL Lacs (yellow curve). FSRQs exhibit a similar behavior up to ${L}_{\gamma }\simeq 3\times {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, where Lν/Lγ strongly increases owing to the increasing neutrino production efficiency because of external radiation fields becoming relevant. It is noteworthy that Lν/Lγ strongly decreases in the range relevant for TXS 0506+056, compared to lower luminosities. This observation will be relevant later on in this analysis.

The right panel of Figure 8 shows the differential contribution to the neutrino flux, where the function dLν/dℓ is obtained as follows:

Equation (10)

Here zmax = 6 for both BL Lacs and FSRQs. Let us recall from Equation (4) that ${10}^{{\ell }}{\epsilon }_{\nu }({\ell })\times \xi ({\ell })={L}_{\nu }({\ell })$, where 10 = Lγ. The factor 1/(1+z)2 refers to the fact that frequency and energy are redshifted. It is clear from this plot that the most important contribution to the high-energy neutrino flux comes from low-luminosity blazars, especially BL Lacs.

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 that the second hump of the SED may be produced by π0 decays. 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 purely leptonic models, since the contribution of these sources to the neutrino flux can be negligible, without affecting the quality of the result, in both scenarios 2 and 3.

4.2. Multiplet Constraint

In this section we discuss the multiplet constraint, referring to the constraint from the nondetection of two or more events from the same source. Namely, we discuss whether the fact that IceCube has not yet observed any multiplet can test our model. This aspect has been emphasized in Murase & Waxman (2016) and Murase et al. (2018), where it was remarked that the absence of neutrino multiplets (or point sources) constrains the contribution of BL Lacs at the level of 10% at 100 TeV and less below 100 TeV. Therefore, in that work blazars are excluded to be the major contributors to the IceCube signal. Our prediction is compatible with the one of Murase et al. (2018) at 100 TeV, since we also find a contribution of ∼10%. On the other hand, we are not interested in explaining the low-energy events, since it is still unclear to which point the background affects events below 100 TeV. We focus instead on the throughgoing muon flux, and therefore on neutrinos above 200 TeV. Under this assumption, the multiplet constraint does not pose a problem to the model, as shown in the remainder of this section.

Let us consider the 9186 sources coming from the distribution by Ajello et al. Each source will contribute with a certain weight to the total neutrino flux expected from blazars, and the weight wi depends on the luminosity Lν and on the redshift z, as follows:

where the relation between Lν and Lγ has been discussed in Sections 3.2 and 3.3 and Dc(z) is defined in Appendix B. Let us remark that more than 90% of the neutrino flux from blazars is powered by BL Lacs with ${L}_{\gamma }\leqslant 3.5\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (about 6300 sources), following our baseline model (scenario 3).

Using the previous information, we can perform several simulations to evaluate the expected number of multiplets, given a certain number of observed signal events. Up to now, 36 throughgoing muons3 have been observed by IceCube (Aartsen et al. 2017c), and 2/3 of them are expected to be signal events, with an uncertainty of about 30%. This translates to a tension between our model and the absence of multiplets of 1σ–1.7σ, as can be seen looking at the orange band of Figure 9. Therefore, the absence of multiplets cannot be considered an issue for our model at this stage.

Figure 9.

Figure 9. Multiplet limit (number of standard deviations σ of exclusion of the model due to the nonobservation of multiplets) as a function of the observed number of signal events in the throughgoing muon data set. In this analysis we use our baseline model (scenario 3) and consider all 9186 blazars contained in the distribution by Ajello et al. We show as an orange band the present number of signals detected using throughgoing muons (Aartsen et al. 2016). The purple band represents the number of years required for IceCube-Gen2 to reach the 5σ level (Aartsen et al. 2014). The band is due to the uncertainty on the contribution of the atmospheric background (mainly prompt neutrinos) to the throughgoing muons, which is about 10% (see Equation (12) of Palladino & Vissani 2017).

Standard image High-resolution image

On the other hand, the next IceCube generation (Aartsen et al. 2014), in which the exposure is expected to be 6–7 times larger than the present one, can test and eventually rule out our model in about 3.5 yr (purple band).

4.3. Comparison with TXS 0506+056

We here discuss the implications of the neutrino associated with TXS 0506+056 in the context of our model.

Expected number of events from TXS 0506+056: The luminosity of neutrinos is expected to be about 1% of the γ-ray luminosity; see middle panel of Figure 8. This is in agreement with the expectations reported in Gao et al. (2018), in which the BL Lac has been modeled using a lepto-hadronic hybrid model. Moreover, we can evaluate the neutrino flux expected from this source, following the same procedure used in Section 2.3 to evaluate the cumulative neutrino flux expected from the entire blazar population. We give three different expected numbers of events, considering the following:

  • 1.  
    The point-source analysis with events above 100 GeV; this data set is dominated by atmospheric background (Aartsen et al. 2016).
  • 2.  
    The throughgoing muon sample with events above 200 TeV; this is dominated by astrophysical signal, with ∼30% of events coming from atmospheric prompt neutrinos (Aartsen et al. 2016).
  • 3.  
    The alert system with extreme high-energy events (EHEs); we will use this to compute the expected number of correlations between neutrinos and known γ-ray sources (Aartsen et al. 2017d).

Using the effective area for point sources reported in Aartsen et al. (2016), we obtain a (time-averaged) expected number of 0.028 events yr−1 from TXS. This is the expected number of events considering the full energy range of the point-source analysis, i.e., with deposited energy larger than 100 GeV. In the left panel of Figure 10 we show the expected number of events per year for objects with the same redshift as TXS but different luminosities. From this plot it is clear that objects with ${10}^{45}\lesssim {L}_{\gamma }(\mathrm{erg}\,{{\rm{s}}}^{-1})\lesssim {10}^{47}$ are the best candidates to produce high-energy neutrinos, due to the fact that for lower luminosities the efficiency epsilon(Lγ) is lower (see Figure 3), whereas for higher luminosities the baryonic loading ξ(Lγ) is too low (see Figure 7) and the neutrino flux is suppressed. To know how many events are expected in the throughgoing muon data set, we must use a low-energy threshold, since only muons above 200 TeV are considered there. A rough estimation can be obtained using Eν = 10 Eμ (see the response function given in Figures S4–S5 of Aartsen et al. 2018b). We therefore define a threshold at 2 PeV for neutrinos, obtaining 0.004 events yr−1 expected in the throughgoing muon data set.

Figure 10.

Figure 10. Left: expected number of events produced by a single source at the same redshift as TXS 0506+056 as a function of the γ-ray luminosity (event rate refers to point-source analysis). The shaded pink region denotes the unresolved sources at redshift z = 0.3365. For a TXS-like source, the (time-averaged luminosity) number is equal to 0.012 events yr−1. Right: distribution of resolved sources that can produce at least the same number of events per year as TXS 0506+056. The different symbols indicate ranges of the expected number of events per year from the source.

Standard image High-resolution image

Let us remark that the observation of one event clearly associated with TXS 0506+056 does not imply that the source has to emit 1 event yr−1. On the contrary, it is possible that there are many faint sources that give a contribution and the detection of one of them is a lucky coincidence. Indeed, the probability of detecting at least one event is given by P(>0) = 1 − exp(−μ), where μ is the expected rate. If μ ≪ 1, the previous expression becomes P(>0) ≃ μ. In the assumption that k sources contribute to the flux, it is sufficient that the rate of a single source is roughly μ ≃ k−1 in order to have a non-negligible probability to observe one event from one source. This concept is known as Eddington bias, and it is explained in Strotjohann et al. (2018) assuming different sources.

Expected number of total events per year: In the right panel of Figure 10, we show the distribution of sources that produce a number of events equal to or greater than that of the TXS source. Using the BL Lac distribution described in Section 2.2, we find that there are 324 visible objects capable of producing at least as many neutrinos per year as TXS. Using the point-source effective area of Aartsen et al. (2018b) and integrating over redshift and luminosity, we find that the average number of events produced by these sources is 0.11 events yr−1 in the point-source analysis (above 100 GeV) and 0.016 events yr−1 in the throughgoing muon data set (above 200 TeV) for each source, which is roughly four times higher than that expected from TXS. Therefore, future correlations will be rather expected at somewhat lower luminosities and redshifts.

In order to estimate how many events per year we expect in IceCube from sources, we have to take into account the following:

  • 1.  
    There are 324 sources capable of producing at least as many neutrinos per year as TXS.
  • 2.  
    Assuming that these 324 sources are isotropically distributed, IceCube can only detect the ones that are visible from the Northern Hemisphere, via throughgoing muons. In particular, the alert system is active for neutrinos above 500 TeV (Aartsen et al. 2017d), and at this energy only neutrinos coming from decl. 0° ≤ δ ≤ 30° can cross the Earth. Therefore, IceCube is roughly sensitive to 1/4 of the sky or, equivalently, to 1/4 of the sources.

Following the previous considerations, the expected number of signal events per year above 100 GeV ${N}_{c}^{100\,\mathrm{GeV}}$ in IceCube is roughly

Equation (11)

These events are hardly distinguishable from the atmospheric background; therefore, it is more interesting to evaluate how many events are expected in the throughgoing muon data set, i.e., above 200 TeV. It is roughly given by

Equation (12)

Let us notice that this number is consistent with the hypothesis that half of the throughgoing muon flux is produced by resolved blazars, since in 8 yr the contribution from these resolved objects would be 10.4 events. At the present, after 8 yr of exposure, 36 throughgoing muons have been observed Aartsen et al. (2017c); 2/3 of them are expected to be signal events, and therefore we expect roughly 12 events from resolved sources.

Expected number of correlations per year: In order to compute the expected number of correlations per year, i.e., neutrinos associated with known objects, we need to refer to the alert system and to use the alert effective area. Using the alert effective area, we obtain 0.004 events yr−1, on average, from each source (including redshift and luminosity distributions). Let us assume a duty cycle of 10%, i.e., the sources are in flaring state only 10% of the time, such as for TXS 0506+056 (Murase et al. 2018). Assuming no neutrino production during the quiescent state, the expected number of events would be 10 times higher during the flaring state, in order to be in agreement with the average result of our work. Therefore, assuming 0.04 events yr−1 during the flaring state and assuming that only 10% of the sources are in the flaring state, the expected number of associations is obtained in the following way:

Equation (13)

The alert system has been active for about 2 yr (Aartsen et al. 2017d), and therefore the expected number of correlations is equal to 0.64, following our baseline model. Up to now, only one correlation has been observed, which is perfectly consistent with our expectation within the Poissonian uncertainty.

4.4. Comparison with PKS 0502+049

For the sake of completeness, we tested the object PKS 0502+049 in the context of our model. In He et al. (2018) the possibility that this object has emitted high-energy neutrinos, during the flaring state observed between 2014 and 2015, is discussed. See also Padovani et al. (2018) for a direct comparison between this object and TXS 0506+056 as neutrino emitters. This source is an FSRQ at redshift z = 0.954, and He et al. (2018) state that to reproduce the γ-ray flare and the neutrino flare the object must exceed the Eddington luminosity by two orders of magnitude during the flaring state. Under this assumption, the flaring luminosity would be ${L}_{\gamma }=9.5\times {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Using the previous luminosity and redshift and under scenario 3, we find that for such an object the neutrino luminosity would be ${L}_{\nu }=7\times {10}^{-5}\,{L}_{\gamma }$ and it would produce a neutrino signal corresponding to the ∼40% of the number of events per year expected from TXS 0506+056. Therefore, the association between the neutrinos observed by IceCube in 2014–2015 and PKS 0502+049 is less plausible than the association between the same neutrinos and TXS 0506+056, but it cannot be ruled out based on our model.

5. Summary and Conclusion

In this work we have studied the possibility that the diffuse astrophysical neutrino flux at the highest energies is powered by blazars, which are AGNs viewed in the direction of the jet. A major obstacle has been the AGN stacking limit, which constrains the blazar contribution from the nonobservation 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 the SED of blazars and γ-ray luminosity, known as the blazar sequence, and a realistic neutrino production model based on these SEDs, we have derived the implications for blazars assuming that the 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, fixing the baryonic loading, high-luminosity objects are very efficient neutrino producers; as a consequence, resolved sources will dominate the neutrino flux, and the stacking limit will constrain 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 a noncontinuous way); in the second one the baryonic loading scales continuously with the γ-ray luminosity and also the expectations from TXS 0506+056 are taken into account. We have found that the only scenario in which blazars can explain the high-energy neutrino flux is that where the baryonic loading of low-luminosity objects is higher than 105, whereas the baryonic loading of high-luminosity sources (both BL Lacs and FSRQs) has to be lower than ∼100. It is also possible that high-luminosity BL Lacs and FSRQs do not contribute to the neutrino flux at all, which means that these sources may not be CR accelerators. Under this hypothesis, low-luminosity objects can then power the entire neutrino flux above a few hundreds of TeV, 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 the constraint coming from the nonobservation of multiplets. While previous works indicate that such high baryonic loadings are difficult to achieve in low-luminosity BL Lacs because of their relatively inefficient accretion, one may speculate that the conditions are different if the neutrinos are produced during flares with baryon injection.

The recent observation of neutrinos from TXS 0506+056 can be interpreted in our baseline scenario. We find that the (time-averaged) expected number of neutrino events from an object at this redshift and luminosity is about 0.004 yr−1 in the throughgoing muon sample. Taking into account the redshift and luminosity distributions, we find 0.016 events per blazar per year on average, coming from about 300 Fermi-LAT-detected objects that are at least as good neutrino emitters as TXS 0506+056. This implies that somewhat larger event rates are expected from objects with higher neutrino luminosity and lower redshifts. Using the effective area of the alert system and taking into account the sky coverage of IceCube, this yields about 0.3 expected neutrino–blazar associations per year. These associations are likely to come from BL Lacs with lower luminosities ${L}_{\gamma }\sim {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and redshifts z ∼ 0.1 in our model.

In conclusion, we have demonstrated that the observed astrophysical flux of throughgoing muons and the stacking limit for blazars are well described only 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. We have not found other possibilities compatible with all the present measurements.

Let us remark that the previous considerations are not necessarily true for a single object, since FSRQs are the most efficient neutrino emitters, as discussed in Rodrigues et al. (2018), and which is evident from Figure 3, meaning that the optical thickness to photohadronic interactions is highest and the protons will transfer energy into neutrinos efficiently. For a diffuse flux, however, this efficiency has to be scaled with the baryonic loading and the number of objects of a certain class. In this work we find that the low-luminosity sources are far more abundant and, given our calculations, must have higher baryonic loadings than the high-luminosity sources if AGNs are to power the diffuse neutrino flux.

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). From the population study, we expect further associations similar to the one to TXS 0506+056 at a rate of about 0.3 yr−1. However, associations with blazars with lower redshift and luminosity than TXS will be more likely.

A.P. thanks F. Vissani for fruitful past discussions concerning the connection between blazars and high-energy neutrinos. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant no. 646623).

Appendix A: Source Model

We present here some details of the source model used in the present work. While most features are similar to 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 2.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 in this work considered to be a soft power law of the blazar luminosity, $B^{\prime} \sim {L}_{\gamma }^{1/5}$. This scaling was employed in Appendix A of Rodrigues et al. (2018) because it can yield values of B' for BL Lacs in the range of 0.1–1 G and for FSRQs in the range of 1 G to a few gauss, 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 B' = 6.3 G for the brightest FSRQ discussed in this work (see Figure 2) and B' = 0.1 G for the dimmest BL Lac. In contrast, in the main text of Rodrigues et al. (2018) a stronger scaling was assumed, $B\sim {L}_{\gamma }^{1/2}$, which yields much weaker magnetic field strengths in the jet of low-luminosity sources.

In Figure 11 we explore the photohadronic 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.

Figure 11.

Figure 11. Left: the solid curves represent the photon density spectra in the jet blob in an FSRQ of luminosity 1048.5 erg s−1 (top) and a BL Lac with 1044.5 erg s−1 (bottom). Note how the external fields of the FSRQ appear boosted in the jet frame owing to the relativistic motion of the jet (see Figure 2). For comparison, we plot in dashed lines the SEDs of an FSRQ and a BL Lac of similar luminosities considered in Rodrigues et al. (2018) (an FSRQ of 1048.8 erg s−1 and a BL Lac of 1044.6 erg s−1). Middle: cooling and interaction rates of protons in the jet for the same sources. The solid red line is the acceleration efficiency considered in this work, η = 10−3 (Equation (1)). For comparison, the red dashed line refers to an ultraefficient proton acceleration (η = 1), discussed in Rodrigues et al. (2018). On the other hand, in that work lower magnetic fields were considered (0.9 G compared to 2 G in this work for the FSRQ and 7 mG compared to 250 mG in this work for the BL Lac). This counteracts the higher acceleration efficiency in the calculation of the acceleration rates and also is reflected in lower synchrotron cooling rates (compare solid and dashed yellow lines). Right: all-flavor neutrino luminosity spectrum produced in the jet, considering the SED and acceleration efficiency used in this work (solid) and in Rodrigues et al. (2018; dashed). Note that in the case of the BL Lac, the suppressed neutrino production below 105 GeV compared to Rodrigues et al. (2018) is simply due to the higher magnetic field, which increases proton synchrotron cooling.

Standard image High-resolution image

In the left panels of Figure 11 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 (${E}_{\gamma }^{{\prime} }={\rm{\Gamma }}{E}_{\gamma }$; see Figure 2). For high-luminosity FSRQs (see end of this section), additional components of the SED are considered in addition to nonthermal radiation (assumed to be produced by electrons in the jet). These external components are an infrared bump from a dusty torus of temperature 3000 K, 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; see Figure 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 Figure 11 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., an acceleration efficiency of η = 10−3 is used (see Equation (1)). In contrast, we show in dashed red lines the case of ultraefficient proton acceleration, as discussed in the main text of Rodrigues et al. (2018). In that work, 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 (see Equation (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 photomeson 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 even if the plasma blob cools more slowly than adiabatically, the size of the blob plays the same role as an adiabatic cooling term in that it limits CR acceleration to the energy where the particle Larmor radius reaches the size of the region, where it cannot be efficiently accelerated further. Note also 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 η = 1.

As we can see from the solid curves in the right panels of Figure 11, the FSRQ (top) produces neutrinos with a maximum of ∼10 PeV in the source frame (∼1 PeV 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 η = 1 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 & Tanaka 2016) indicate that high acceleration efficiencies may be difficult to achieve in relativistic shocks, as is the case for blazar jets. On the other hand, with such low acceleration efficiencies, these sources cannot power the diffuse ultra-high-energy CR spectrum (see bottom panel of Figure 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 1-day variability timescale, ${r}_{\mathrm{blob}}^{{\prime} }={\rm{\Gamma }}c\times 1\,\mathrm{day}=3\times {10}^{16}\,\mathrm{cm}$, other values could be considered. In low-luminosity blazars, the maximum neutrino energy then scales as ${E}_{\nu }^{\max }\sim {r}_{\mathrm{blob}}^{{\prime} }$, because CR acceleration is limited only by the blob size (see bottom middle panel of Figure 11); in high-luminosity objects, this scaling is weaker because the maximum energy is limited by photohadronic interactions (see top middle panel of Figure 11). The total neutrino luminosity produced by the source scales mainly with the optical thickness to photomeson production. The optical thickness grows linearly with ${r}_{\mathrm{blob}}^{{\prime} }$ and with ${n}_{\gamma }\sim V{{\prime} }^{-1}\sim r{{\prime} }_{\mathrm{blob}}^{-3};$ thus, the normalization of the emitted neutrino spectrum scales as $r{{\prime} }_{\mathrm{blob}}^{-2}$.

Finally, the model geometry (briefly discussed in Section 2.1) follows the principles discussed in detail in Rodrigues et al. (2018). The jet blob is at a distance to the black hole ${r}_{\mathrm{diss}}={r}_{\mathrm{blob}}^{{\prime} }/\sin ({\theta }_{\mathrm{jet}})={\rm{\Gamma }}{r}_{\mathrm{blob}}^{{\prime} }$ (see Figure 1), which corresponds to the dissipation radius of the jet, assumed the same for all blazars. At the same time, the sizes of the BLR and dusty torus (rBLR and rDT) are assumed to scale with ${L}_{\mathrm{disk}}^{1/2}$ (Murase et al. 2014; Rodrigues et al. 2018), where Ldisk is the accretion disk luminosity. In turn, Ldisk and Lγ also scale together as ${L}_{\gamma }\propto {L}_{\mathrm{disk}}^{0.683}$ (Murase et al. 2014; Rodrigues et al. 2018). Therefore, the jet blob lies inside the BLR only in FSRQs with luminosity ${L}_{\gamma }\gt 2.4\times {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. 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 rBLR < rdiss < rDT, i.e., with luminosity $3\times {10}^{46}\lt {L}_{\gamma }(\mathrm{erg}\,{{\rm{s}}}^{-1})\lt 2.4\times {10}^{48}$, 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 Figure 2 in which the thermal IR bump is present, but no optical bumps or emission lines). This approximation is justified by the conclusion in Rodrigues et al. (2018) that for a purely diffusive CR escape, the neutrinos produced during CR propagation through the BLR typically contribute negligibly to the total neutrino fluence. Rather, the flux is dominated either by the neutrinos produced directly in the jet or during CR propagation through the thermal field from the dust torus.

Appendix B: The Distribution of BL Lacs and FSRQs

The parameterization in Ajello et al. (2014, 2012) describes the differential source distribution

Equation (14)

which represents the number of sources N per comoving volume Vc, Lγ is the emitted luminosity between 0.1 and 100 GeV, and Γ is the slope of the γ-ray flux, assuming power-law spectra in that energy interval.

For our purposes, it is convenient to write the previous expression as a function of the redshift z and ${\ell }\equiv {\mathrm{log}}_{10}[\tfrac{{L}_{\gamma }}{\mathrm{erg}\,{{\rm{s}}}^{-1}}]$. The comoving volume is given by

Equation (15)

The comoving distance Dc(z) is defined as follows:

Equation (16)

where H0 = 4.22 Gpc is the Hubble distance. The function d(z) is given by

Equation (17)

where $h(z)=\sqrt{{{\rm{\Omega }}}_{{\rm{\Lambda }}}+{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}}$, with ΩΛ = 0.73 and Ωm = 0.27. Based on the previous equations, the relation between the comoving volume and the redshift is

Equation (18)

The Jacobian obtained from the transformation ${L}_{\gamma }\to {\ell }$ is given by

Equation (19)

Combining the two previous expressions with Equation (14), we obtain

Equation (20)

The parameterization by Ajello et al. is characterized by several parameters, which we report in Table 2. Only the first two parameters, namely, A and L*, are dimensional, in Gpc−3 and erg s−1, respectively. In this discussion, the quantity Lγ is always in units of erg s−1.

Table 2.  Table of the Parameters That Define the Distribution of BL Lacs and FSRQs

  A(Gpc−3) L* (erg s−1) γ1 γ2 μ* β σ ${p}_{1}^{* }$ ${p}_{2}^{* }$ τ ${z}_{c}^{* }$ α
BL Lacs 3.39 1047.4472 0.27 1.86 2.10 0.0646 0.26 2.24 −7.37 4.92 1.34 0.0453
FSRQs 3.06 1047.9243 0.21 1.58 2.44 0 0.18 7.35 −6.51 0 1.47 0.21

Note. Only the first two parameters are dimensional, as indicated.

Download table as:  ASCIITypeset image

The distribution in the slope Γ is assumed to be Gaussian,

Equation (21)

where $\mu ({L}_{\gamma })={\mu }_{* }+\beta [{\mathrm{log}}_{10}({L}_{\gamma })-46]$.

The distribution in luminosity is given by a double power law:

Equation (22)

The distribution in z is described by an evolution factor, also a double power law:

Equation (23)

where ${z}_{c}{({L}_{\gamma })={z}_{c}^{* }\times ({L}_{\gamma }/{10}^{48})}^{\alpha }$ and ${p}_{1}({L}_{\gamma })={p}_{1}^{* }\,+\tau [{\mathrm{log}}_{10}({L}_{\gamma }-46)]$. Note: in Ajello et al. (2012, 2014) the previous equation is reported with a wrong positive sign for both p1 and p2. This mistake has been corrected later, in a footnote in Ajello et al. (2015). The wrong sign produces an important difference, especially in the distribution of low-luminosity BL Lacs, which accumulates at high redshift using the wrong distribution.

Summarizing, we have

Equation (24)

which, using Equation (20), becomes

Equation (25)

The coefficient ${ \mathcal N }=4\pi {D}_{H}^{3}\times A$ is a dimensionless number.

Our goal is to compute the diffuse high-energy neutrino flux from blazars. Since the neutrino spectra at the source only depend on the luminosity, we can integrate over the slope Γ (we cannot integrate over the redshift z, 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:

Equation (26)

The intervals of integration for BL Lacs are z ∈ [0.03,6] and  ∈ [43.85,52], and those for FSRQs are z ∈ [0,6] and  ∈ [44,52]. The distributions $\tfrac{{d}^{2}N}{{dzd}{\ell }}$ of BL Lacs and FSRQs are represented in the left panel of Figure 4.

Footnotes

  • Primed 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 Ajello et al. (2014, 2012) Equation (23) is reported with the wrong signs of p1 and p2. This produces an incorrect distribution, which is particularly relevant for low-luminosity BL Lacs. In this work we use this function considering the correct signs. Please note that the erratum has been reported in a later work in 2015. See the second footnote of Ajello et al. (2015).

  • In order to evaluate the expected number of multiplets, we need events with good angular resolution; therefore, we use the throughgoing muon data set in this analysis since the average angular resolution is ∼1°.

Please wait… references are loading.
10.3847/1538-4357/aaf507