Proton Synchrotron Gamma-Rays and the Energy Crisis in Blazars

High-energy emission of blazars—active galactic nuclei (AGN) with relativistic jets closely aligned to our line of sight, powered by accretion onto a supermassive black hole (BH)—has been a matter of vibrant debate since their first detection in γ-rays (for a review, see Blandford et al. 2018).

Historically, γ-ray emission has been attributed to two broad classes of models that are distinguished mainly by the species of radiating particles. Leptonic models invoke inverse Compton scattering of low-energy photons by relativistic electrons (e.g., Marscher & Gear 1985; Dermer et al. 1992). Hadronic models involve a variety of mechanisms that are directly or indirectly related to relativistic hadrons, such as proton synchrotron (PS) radiation (e.g., Aharonian 2000), or synchrotron and Compton processes of secondary electrons and positrons produced in photohadronic interactions (e.g., Mannheim 1993).

Unveiling the dominant process for blazar's γ-ray emission has been the subject of numerous studies. This is not surprising, as by constraining the dominant high-energy processes in blazars we can probe the jet's physical conditions (which are hidden to direct observation) and help answer longstanding questions regarding launching and mass-loading of jets (Blandford et al. 2018).

The most common methods to probe the origin of γ-rays are spectral energy distribution (SED) modeling of broadband emission (e.g., Böttcher et al. 2013; Ghisellini & Tavecchio 2015; Petropoulou et al. 2015) and searches for correlated variability between low-energy radiation (e.g., radio and optical) and γ-rays (e.g., Max-Moerbeck et al. 2014; Liodakis et al. 2018b, 2019). While past studies have favored leptonic models, they have not been always conclusive. The most recent possible association of high-energy neutrinos with blazar TXS 0506+056 (IceCube Collaboration et al. 2018a, 2018b) would also suggest that the usually disfavored hadronic component should be present. Interestingly, SED modeling of the first likely multi-messenger event point to leptonic processes dominating the γ-ray emission (e.g., Gao et al. 2019; Zhang et al. 2020), although the jet's energetics are still governed by relativistic hadrons (Keivani et al. 2018; Petropoulou et al. 2020).

Indeed, one of the major criticisms of hadronic models for blazar emission relates to their energetic requirements. The inefficiency of hadronic processes was pointed out using generic arguments by Sikora et al. (2009), Sikora (2011), and later discussed on a source-to-source basis using SED modeling of steady emission or γ-ray flares (e.g., Cerruti et al. 2015; Petropoulou & Dimitrakoudis 2015; Petropoulou et al. 2017). Recently, Zdziarski & Böttcher (2015) explored the energetic requirements of the PS model for a limited sample (12 sources from Böttcher et al. 2013), and concluded that the estimated minimum jet powers are not compatible with the inferred accretion power and Eddington luminosity.

In light of recent results, we revisit the jet-power analysis of γ-ray blazars in the PS model by following the analytical approach of Petropoulou & Dermer (2016) and extending our calculations to the largest sample to date (145 sources).

Our sample consists of sources with synchrotron peak frequency and luminosity from the fourth Fermi AGN catalog (4LAC; The Fermi-LAT Collaboration 2019), Doppler factors from Liodakis et al. (2018a), and apparent velocities (β_app) from the MOJAVE survey Lister et al. (2016). We use the SED builder tool³ of the Space Science Data Center (SSDC) to fit a third-degree polynomial in log–log space to the archival data (similar to the analysis of the 4LAC for the synchrotron spectrum) to estimate the peak frequency and luminosity of the high-energy component. We have removed sources where the data are insufficient to confidently determine the parameters of the high-energy component; these sources either lacked X-ray observations that constrain the low-energy part of the high-energy component or γ-ray observations were not available through SSDC. Our final sample consists of 145 sources.

The absolute jet power for a two-sided jet can be written as ${P}_{j}=2\pi {r}^{{\prime} 2}\beta {{\rm{\Gamma }}}^{2}c{\sum }_{i=B,e,p}\left({u}_{i}^{{\prime} }+{p}_{i}^{{\prime} }\right)+{P}_{j}^{r}+{P}_{j}^{c}$, where r' is the radius of the emitting region⁴ , Γ = (1 − β²)^−1/2 is the jet's bulk Lorentz factor, ${u}_{i}^{{\prime} }$ is the energy density of relativistic particles/magnetic fields and ${p}_{i}^{{\prime} }={u}_{i}^{{\prime} }/3$, P_j^r is the absolute photon luminosity, and P_j^c is the contribution of cold protons to the total jet power (Zdziarski 2014; Petropoulou & Dermer 2016). Henceforth, we drop the latter term from our analysis, as it is negligible compared to the others in the PS scenario.

Following Petropoulou & Dermer (2016, hereafter PD16), we assume monoenergetic particle distributions for both relativistic electrons and protons (i.e., ${N}_{i}({\gamma }_{i}^{{\prime} })\,={N}_{i,0}\delta ({\gamma }_{i}^{{\prime} }-{\bar{\gamma }}_{i}^{{\prime} })$, $i=e,p$). This choice is equivalent to the assumption of power-law particle energy spectra with slopes p < 2, while consideration of steep power laws (p > 2) would only increase our minimum power estimates. As in PD16, we assume that the proton radiative efficiency is ≃1 (lower efficiency would only increase the energetic requirements). We can re-write P_j as a function of the emitting region's Doppler factor ${ \mathcal D }={\left[{\rm{\Gamma }}(1-\beta \cos \theta )\right]}^{-1}$ (here, θ is the observer's angle) and co-moving magnetic field strength x ≡ B'/B_cr (in units of B_cr = 4.4 × 10¹³ G),

$$\begin{eqnarray}&&\begin{array}{l}{\psi }^{-2}{P}_{j}={A}_{B}({t}_{{\rm{v}}},z){x}^{2}{{ \mathcal D }}^{4}\\ \quad +\,{A}_{e}({L}_{l},{\varepsilon }_{l},{t}_{{\rm{v}}},z){x}^{-3/2}{{ \mathcal D }}^{-5/2}\left(1+\displaystyle \frac{2{A}_{e}}{{L}_{l}}{x}^{-3/2}{{ \mathcal D }}^{-1/2}\right)\\ \quad +\,{A}_{p}({L}_{h},{\varepsilon }_{h},{t}_{{\rm{v}}},z){x}^{-3}{{ \mathcal D }}^{-4}\left(1+\displaystyle \frac{2{A}_{p}}{{L}_{h}}{x}^{-3}{{ \mathcal D }}^{-2}\right)\\ \quad +\,{A}_{r}({L}_{h}/{L}_{l}){{ \mathcal D }}^{-2}{L}_{l}.\end{array}\end{eqnarray} \tag{ 1 }$$

Here, $\psi =1+{({\rm{\Gamma }}\theta )}^{2}\approx 2{\rm{\Gamma }}/{ \mathcal D }$, and ${A}_{i}(\cdots )$ with $i=B,e,p,r$ are functions of source parameters⁵ : redshift z, typical variability timescale ${t}_{{\rm{v}}}$, peak luminosities of the low- and high-energy SED humps, L_l and L_h, and the respective peak photon energies ε_l and ε_h (both in units of m_e c²). Knowledge of the SED parameters, relativistic boosting effects, and variability timescale of a source allows us to estimate the minimum jet power (${P}_{j,\min }$) with respect to the unknown variable B'.

For each source we derive ${P}_{j,\min }$ and the corresponding magnetic field strength ${B}^{{\prime} }$ for 10⁴ combinations of random values for ${\varepsilon }_{l,h}$, ${L}_{l,h}$, ${t}_{{\rm{v}}}$, ${ \mathcal D }$, and ${\beta }_{\mathrm{app}}$ drawn from Gaussian distributions with mean μ and standard deviation σ. For ${t}_{{\rm{v}}}$ we choose μ = 10⁵ s and σ = 3 × 10⁴ s, which translates to a range of minutes to > 2 days (e.g., Meyer et al. 2019). We assume a σ of 0.5 dex for the luminosities and for the peak frequencies 0.3 dex (Lister et al. 2015). For ${ \mathcal D }$ and β_app needed to estimate Γ and θ, we use the values and their uncertainties listed in Liodakis et al. (2018a).

To assess our results, we compare the derived ${P}_{j,\min }$ to three characteristic "energy estimators" of an accreting BH system: (i) the Eddington luminosity ${L}_{\mathrm{Edd}}$, (ii) the accretion disk luminosity L_d, and (iii) the power of the Blandford–Znajek (BZ) process P_BZ (Blandford & Znajek 1977). We estimate the BH masses for 82 blazars in our sample (needed for computing L_Edd and P_BZ) using the Hβ, Mg ii, and C iv FWHM and line luminosities from Shaw et al. (2012) and Torrealba et al. (2012), together with the scaling relations from Shaw et al. (2012, Equation (5), Table 2). We complement our sample with 13 mass estimates from Woo & Urry (2002), Wang et al. (2004), and Liu et al. (2006) that use the same lines. For the remaining sources we use the BL Lac and Flat Spectrum Radio Quasar (FSRQ) population median and standard deviations derived from the BH estimates in this work. The accretion disk luminosity (71 sources) is estimated using the line luminosities of the Hβ, Mg ii, and O iii lines and the scaling relations from Zamaninasab et al. (2014, Equations (9)–(11)). When multiple estimates for either the BH mass or L_d are available we use the median of the estimates and for its uncertainty we quote the standard deviation of the estimates or the average uncertainty (whichever is greater).

To estimate the power of the BZ process, we first estimate the jet's co-moving magnetic field strength at 1 parsec (pc), ${B}_{1\mathrm{pc}}^{{\prime} }$, using the core-shift measurements and Equations (2) and (3) from Zamaninasab et al. (2014), with the correct redshift terms (Lobanov 1998; Zdziarski et al. 2015). We then derive the poloidal magnetic flux that threads the pc-scale jet, Φ_jet, using the jet apparent opening angles from Pushkarev et al. (2009) and Equation (1) from Zamaninasab et al. (2014). This quantity is a proxy of the poloidal magnetic flux threading the BH (Φ_BH) under the flux-freezing assumption. For the BH spin a, we consider three cases: all BHs are maximally spinning; a follows a uniform distribution from 0 to 1; a follows a Beta distribution⁶ with μ = 0.937, σ = 0.074 for BL Lacs and μ = 0.742, σ = 0.163 for FSRQs (Liodakis 2018). The BZ power is then estimated as ${P}_{\mathrm{BZ}}=\kappa {{\rm{\Omega }}}_{{\rm{H}}}^{2}{{\rm{\Phi }}}_{\mathrm{BH}}^{2}f({{\rm{\Omega }}}_{{\rm{H}}})/4\pi c$, where $\kappa \approx 0.05$ is a numerical constant whose value depends on the magnetic field geometry, ${{\rm{\Omega }}}_{{\rm{H}}}={ac}/2{r}_{{\rm{H}}}$ is the angular frequency of the BH horizon, ${r}_{{\rm{H}}}={r}_{g}(1+\sqrt{1-a})$ is the BH event horizon radius, ${r}_{g}={{GM}}_{\bullet }/{c}^{2}$ is the gravitational radius, c is the speed of light, M_• is the BH mass, and $f({{\rm{\Omega }}}_{{\rm{H}}})$ as $f{({{\rm{\Omega }}}_{{\rm{H}}})\approx 1+1.38({{\rm{\Omega }}}_{{\rm{H}}}{r}_{g}/c)}^{2}-9.2{({{\rm{\Omega }}}_{{\rm{H}}}{r}_{g}/c)}^{4}$ (Tchekhovskoy et al. 2011). All derived parameters are listed in Table 1.

Table 1.  Parameter Estimates for the Sources in Our Sample

Name	Alt. Name	${\epsilon }_{h}$	Lh	${B}^{{\prime} }$	${P}_{j,\min }$	M•	Ld	PBZ	${\epsilon }_{\nu }^{{pk}}$	${\epsilon }_{\nu }^{{pk}}{F}_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}^{{pk}}$
4FGLJ0017.5-0514	J0017-0512	0.0332	45.05	16.88${}_{-5.77}^{+14.75}$	47.14${}_{-0.21}^{+0.38}$	7.55 ± 0.45	45.92 ± 0.51	45.86 ± 1.24	5.07${}_{-0.19}^{+0.29}$	−15.15${}_{-0.34}^{+1.64}$
4FGLJ0019.6+7327	J0019+7327	0.0605	48.02	187.38${}_{-64.10}^{+97.42}$	48.42${}_{-0.20}^{+0.38}$	9.62 ± 0.52	47.15 ± 0.52	46.15 ± 1.21	4.43${}_{-0.16}^{+0.24}$	−13.96${}_{-0.32}^{+1.38}$
4FGLJ0051.1-0648	J0051-0650	0.0937	47.27	168.76${}_{-95.71}^{+311.88}$	48.82${}_{-0.38}^{+4.31}$	⋯	⋯	46.86 ± 1.48	4.47${}_{-0.26}^{+0.46}$	−13.55${}_{-0.38}^{+5.07}$
4FGLJ0108.6+0134	J0108+0135	0.1039	47.88	81.11${}_{-46.00}^{+175.39}$	49.26${}_{-0.30}^{+2.39}$	⋯	⋯	47.24 ± 1.36	4.77${}_{-0.26}^{+0.56}$	−13.66${}_{-0.39}^{+4.84}$
4FGLJ0112.8+3208	J0112+3208	0.0147	46.44	59.26${}_{-27.63}^{+64.03}$	48.99${}_{-0.24}^{+1.02}$	⋯	⋯	46.94 ± 1.33	4.55${}_{-0.23}^{+0.38}$	−14.02${}_{-0.36}^{+2.74}$
4FGLJ0116.0-1136	J0116-1136	0.0093	46.02	73.05${}_{-34.07}^{+95.71}$	48.94${}_{-0.26}^{+1.09}$	8.77 ± 0.38	45.32 ± 1.08	46.90 ± 1.36	4.35${}_{-0.23}^{+0.39}$	−14.34${}_{-0.36}^{+2.74}$
4FGLJ0132.7-1654	J0132-1654	0.0178	46.75	53.37${}_{-24.89}^{+57.67}$	49.34${}_{-0.27}^{+0.97}$	⋯	⋯	47.36 ± 1.29	4.61${}_{-0.21}^{+0.37}$	−14.28${}_{-0.35}^{+2.34}$
4FGLJ0137.0+4751	J0136+4751	0.0937	46.86	15.20${}_{-9.27}^{+57.85}$	48.40${}_{-0.21}^{+0.77}$	8.68 ± 0.31	46.23 ± 0.52	46.17 ± 1.07$\dagger$	5.35${}_{-0.31}^{+0.63}$	−14.56${}_{-0.39}^{+10.68}$
4FGLJ0152.2+2206	J0152+2207	0.0314	46.49	187.38${}_{-87.38}^{+293.26}$	49.12${}_{-0.30}^{+2.21}$	⋯	⋯	46.98 ± 1.39	4.19${}_{-0.24}^{+0.42}$	−13.31${}_{-0.37}^{+3.54}$
4FGLJ0204.8+1513	J0204+1514	0.0161	46.27	168.76${}_{-57.73}^{+116.04}$	49.36${}_{-0.24}^{+0.60}$	⋯	⋯	47.21 ± 1.24	4.15${}_{-0.19}^{+0.28}$	−13.92${}_{-0.33}^{+1.53}$

Note. The values of all parameters except ${\epsilon }_{h}$ and ${B}^{{\prime} }$ are displayed as log₁₀. ${\epsilon }_{h}$ is in GeV, ${\epsilon }_{\nu }^{{pk}}$ in TeV, ${B}^{{\prime} }$ in Gauss, M_• in solar masses, ${\epsilon }_{\nu }^{{pk}}{F}_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}^{{pk}}$ in TeV cm⁻² s⁻¹, and L_h, ${P}_{j,\min }$, L_d, ${{\rm{P}}}_{\mathrm{BZ}}$ in erg s⁻¹. All ${P}_{\mathrm{BZ}}$ estimates derived from core-shift measurements are indicated with $\dagger$.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Panels (a) and (b) of Figure 1 show, respectively, the distributions of ${P}_{j,\min }$ and ${B}^{{\prime} }$ that minimizes the total jet power for 3C 273. The magnetic field and relativistic proton components contribute the most to the total jet power, as expected in the PS scenario (PD16). Our analytical method yields ${P}_{j,\min }={8.6}_{-4.6}^{+11.7}\times {10}^{49}$ erg s⁻¹ which is consistent with SED modeling results (Böttcher et al. 2013; Petropoulou & Dimitrakoudis 2015). We also find comparable (within 90% uncertainty) ${P}_{j,\min }$ values for 10 sources we have in common with Böttcher et al. (2013). Panel (c) of Figure 1 shows the results for the whole sample. Different symbols are used to identify blazar classes according to their peak (rest-frame) synchrotron frequency: low-synchrotron peaked (LSP) sources (${\nu }_{s}\lt {10}^{14}$ Hz), intermediate-synchrotron peaked (ISP; ${10}^{14}\,\lt {\nu }_{s}\lt {10}^{15}$ Hz), and high-synchrotron peaked (HSP; ν_s > 10¹⁵ Hz) sources (Abdo et al. 2010). The minimum jet power decreases on average as we move from LSP to HSP sources (PD16), while blazars with higher ${P}_{j,\min }$ tend to have stronger magnetic fields in the emission region.

Zoom In Zoom Out Reset image size

Figure 2 shows the comparison of the minimum jet power with ${L}_{\mathrm{Edd}}$ (top panel) and L_d (bottom panel). None of the ISP and HSP sources in our sample have BH masses, thus their L_Edd is computed using the population estimates (Section 3). Except for a handful of sources with ${P}_{j,\min }\sim {L}_{\mathrm{Edd}}$ (within uncertainties), we find that the majority of blazars in the PS scenario has super-Eddington jet powers and ${P}_{j,\min }\sim {10}^{2}{L}_{d}$ (see also Zdziarski & Böttcher 2015).

Zoom In Zoom Out Reset image size

Figure 3 (top panel) shows the comparison of ${P}_{j,\min }$ with ${P}_{\mathrm{BZ}}$ for 40 blazars with core-shift measurements (open colored symbols), assuming that all sources host maximally spinning BHs. We have also estimated the BZ power for sources without core-shift measurements (filled gray symbols) using the sample's median (and standard deviation) opening angle and magnetic field. Most sources cluster around the ${P}_{j,\min }={10}^{2}{P}_{\mathrm{BZ}}$ line, and the deviation from the line of equality becomes even larger when considering uniform or beta distributions for the BH spin (not shown in the figure). Meanwhile, we find that ${{\rm{\Phi }}}_{\mathrm{BH}}\approx 50{\left(\dot{M}{r}_{g}^{2}c\right)}^{1/2}\propto {L}_{d}^{1/2}{M}_{\odot }$ (bottom panel), as expected for magnetically arrested accretion disks (Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan et al. 2003), in agreement with Zamaninasab et al. (2014). Thus, the PS scenario predicts much higher jet powers than the BZ power, even in the MAD regime where the jet production efficiency is highest (Tchekhovskoy et al. 2011). While measured ${P}_{\mathrm{BZ}}$ are only available for 40 sources in our sample, the on-average estimates of the remaining sources follow the same trend well. Equation (1) from Zamaninasab et al. (2014) used to estimate Φ_jet, assumes energy equipartition between magnetic fields and radiating particles and does not explicitly consider the relation between the jet opening angle and magnetization σ_M. By relaxing this assumption, Zdziarski et al. (2015) derived a more general expression (see their Equation (21)), which is identical to that of Zamaninasab et al. (2014) for σ_M = 1, but yields lower Φ_jet values (by a factor of 2^−1/2) for σ_M ≪ 1. While a small correction given the uncertainty of individual estimates, it would only increase the discrepancy between P_BZ and P_j,min.

Zoom In Zoom Out Reset image size

Because of several assumptions made in this work (i.e., Doppler factor estimates, monoenergetic particle distribution, and proton radiative efficiency), the derived values of ${P}_{j,\min }$ constitute lower limits of the true minimum jet power further increasing this discrepancy. Hence, our results strongly disfavor the PS scenario for the majority of blazars, particularly for LSPs.

Location of γ-ray emission region. We can estimate the location of the γ-ray production site for the sources having estimates of the pc-scale jet's magnetic field as follows. Assuming that the jet's magnetic field is roughly equal to the magnetic field strength of the emission region, i.e., ${B}^{{\prime} }\approx {B}_{j,\phi }^{{\prime} }\propto 1/z$, we may write ${z}_{\mathrm{em}}\approx ({B}_{1\mathrm{pc}}^{{\prime} }/{B}^{{\prime} })\,\mathrm{pc}$. We then find that ${z}_{\mathrm{em}}\ll 1$ pc, with 68% of the values ranging between 0.006 and 0.08 pc, with a median of 0.03 pc. Given that the median radius of the broad line region (BLR) for the sources of our sample is 0.15 pc, our results suggest that the γ-ray production site should be well within the BLR. This conclusion is, however, in tension with the lack of strong absorption features in the GeV γ-ray spectrum of luminous quasars (e.g., Costamante et al. 2018). The sub-pc location of the emission region is also inconsistent with the radius inferred by the average observed variability, i.e., ${r}^{{\prime} }=c{ \mathcal D }{t}_{{\rm{v}}}/(1+z)$. The cross-sectional radius of the jet at the emission region can be written as ${\varpi }_{\mathrm{em}}\approx {z}_{\mathrm{em}}{\theta }_{j}$, for a conical jet with small half-opening angle θ_j (the same assumption is made when computing ${B}_{1\mathrm{pc}}^{{\prime} }$). Although a consistent picture would require ${r}^{{\prime} }\lesssim {\varpi }_{\mathrm{em}}$, we find ${r}^{{\prime} }/{\varpi }_{\mathrm{em}}\gt 1$, with 68% of the ratio values in the range 9–60 and a median of 27.

Part of this tension can be resolved, if one assumes that the magnetic field in the emission region is amplified with respect to the jet's toroidal magnetic field component. By writing ${B}^{{\prime} }={f}_{\mathrm{amp}}{B}_{j,\phi }^{{\prime} }$ and requiring ${r}^{{\prime} }={\varpi }_{\mathrm{em}}$, we find that the median amplification factor needed is 27. Thus, the γ-ray production site is also moved to pc scales, typically beyond the BLR (median z_em = 0.5 pc and 68% of values ranging between 0.2 and 1.6 pc). Alternatively, lower B' values can be derived if the emission region moves with larger ${ \mathcal D }$ than what we have assumed (e.g., a three times higher ${ \mathcal D }$ for all sources would yield ${B}^{{\prime} }\sim 1\mbox{--}100$ G), but at the cost of higher ${P}_{j,\min }$.

High-energy neutrino emission. Relativistic protons can also interact with low-energy photons to produce high-energy electron and muon neutrinos through the photomeson ($p\pi$) production process. The apparent isotropic proton luminosity L_p and the absolute jet power in relativistic protons ${P}_{j,p}$ are related as ${L}_{p}=2{{ \mathcal D }}^{2}{\psi }^{-2}{P}_{j,p}$ (e.g., Dermer et al. 2012). For the purposes of this discussion, we replace the mononenergetic proton distribution with a power-law spectrum with an exponential cutoff, so that the differential apparent isotropic proton luminosity is written as ${L}_{p}({\epsilon }_{p})\propto {\epsilon }_{p}^{-p}{e}^{-{\epsilon }_{p}/{\epsilon }_{p,\max }}$, where p = 1.7 and ${\epsilon }_{p,\max }={ \mathcal D }{\bar{\gamma }}_{p}^{{\prime} }{m}_{p}{c}^{2}/(1+z)$. For every source, we compute ${L}_{p}({\epsilon }_{p})$ and ${\epsilon }_{p,\max }$ for parameters minimizing the total jet power (Sections 3 and 4).

Following our previous discussion on the location of the emission region, we assume that protons interact only with the jet's synchrotron photons. The differential number density of the low-energy photons is ${n}^{{\prime} }(x)={n}_{0}^{{\prime} }{\varepsilon }_{l}^{{\prime} -2}\left[{x}^{-2+{{\rm{\Gamma }}}_{1}}H(1-x)+{x}^{-2+{{\rm{\Gamma }}}_{2}}H(x-1)\right]$, where ${{\rm{\Gamma }}}_{1}=1/2$, ${{\rm{\Gamma }}}_{2}\,=-1/2$ ⁷ , $x\equiv {\varepsilon }^{{\prime} }/{\varepsilon }_{l}^{{\prime} }$, ${\varepsilon }_{l}^{{\prime} }={\varepsilon }_{l}(1+z)/{ \mathcal D }$, and ${n}_{0}^{{\prime} }=3{L}_{l}{(1+z)}^{2}/4\pi {m}_{e}{c}^{5}{t}_{{\rm{v}}}^{2}{{ \mathcal D }}^{6}$. The $p\pi$ efficiency is defined as ${f}_{p\pi }\equiv { \mathcal D }{t}_{{\rm{v}}}/{t}_{p\pi }^{{\prime} }(1+z)$, where ${t}_{p\pi }^{{\prime} }$ is the energy-loss timescale. This is ${t}_{p\pi }^{{\prime} -1}\left({\gamma }_{p}^{{\prime} }\right)\,=\,c/(2{\gamma }_{p}^{{\prime} 2}){\int }_{0}^{\infty }d\varepsilon ^{\prime} n^{\prime} (\varepsilon ^{\prime} )/{\varepsilon }^{{\prime} 2}{\int }_{{\varepsilon }_{\mathrm{th}}}^{2{\gamma }_{p}^{{\prime} }{\varepsilon }^{{\prime} }}d{\varepsilon }_{r}{\sigma }_{p\pi }({\varepsilon }_{r}){\kappa }_{p\pi }({\varepsilon }_{{\rm{r}}}){\varepsilon }_{r}$ (Stecker 1968), where ${\varepsilon }_{\mathrm{th}}\approx 400$ is the threshold photon energy for production of a ${{\rm{\Delta }}}^{+}(1232)$ resonance, ${\kappa }_{p\pi }=0.2$ is the inelasticity of interaction, and ${\sigma }_{p\pi }\approx 0.34\,\mathrm{mb}$ for ${\varepsilon }_{\mathrm{th}}\leqslant {\varepsilon }_{r}\,\leqslant 980$ is the cross section (Dermer & Menon 2009). The differential all-flavor neutrino (and anti-neutrino) flux is given by

$$\begin{eqnarray}&&{\epsilon }_{\nu }{F}_{\nu +\bar{\nu }}({\epsilon }_{\nu })\approx \displaystyle \frac{3}{8}{f}_{p\pi }\left(\displaystyle \frac{{\epsilon }_{p}(1+z)}{{ \mathcal D }{m}_{p}{c}^{2}}\right)\displaystyle \frac{{\epsilon }_{p}{L}_{p}({\epsilon }_{p})}{4\pi {d}_{L}^{2}},\end{eqnarray} \tag{ 2 }$$

where ${\epsilon }_{\nu }\approx {\epsilon }_{p}/20$ and d_L is the luminosity distance. Because of neutrino oscillations the muon neutrino and anti-neutrino energy flux at Earth is ${F}_{{\nu }_{\mu }+{\bar{\nu }}_{\mu }}\approx {F}_{\nu +\bar{\nu }}/3$.

Figure 4 shows the peak neutrino energy and peak ${\nu }_{\mu }+{\bar{\nu }}_{\mu }$ energy flux derived from Equation (2). Our results are in line with predictions made for individual sources, i.e., ∼0.1–1 EeV neutrinos with fluxes much lower than in γ-rays (e.g., Dimitrakoudis et al. 2014; Keivani et al. 2018). There are a few blazars that are potentially interesting neutrino sources (close to IceCube's discovery potential), with LSP blazar 4FGLJ2148.6+0652 being the best example.

Zoom In Zoom Out Reset image size

For this blazar, we find ${B}^{{\prime} }\sim 2565$ G and ${P}_{j,\min }\gtrsim {10}^{2}{L}_{d}$, which is roughly three orders of magnitude higher than the average P_BZ value of LSPs with core-shift measurement. If steady neutrino emission at the predicted flux levels is detected from this blazar by IceCube, or from other sources by future experiments, such as GRAND (Álvarez-Muñiz et al. 2020) or POEMMA (Venters et al. 2019) our understanding of accretion and jet launching in blazars needs to be revised.

We explored the energetic requirements for the PS model for the largest sample of γ-ray blazars. The expectation for the minimum jet power in our sample far exceeds the observed L_d and L_Edd as well as the derived P_BZ, even more so, when considering that the results of this work constitute a lower limit to the true minimum jet power. In addition, the derived magnetic field strengths in the emission region imply either large amplification of the jet's magnetic field or sub-pc γ-ray production sites, well within the BLR and in tension with recent results. The expected neutrino emission (for all sources in our sample) peaks at ∼0.1–10 EeV, i.e., at much higher energies than the multi-TeV neutrinos associated with TXS 0506+056 (see also Keivani et al. 2018). Meanwhile, the typical peak neutrino fluxes are ∼10⁻⁴ times lower than the peak γ-ray fluxes. Our results clearly demonstrate that the scenario where PS accounts for the observed steady γ-ray emission in blazars is highly unlikely. Given that alternative hadronic models invoking emission from $p\pi$ secondaries typically require even higher-energy budgets (e.g., Petropoulou et al. 2015), we conclude that if a hadronic population is present in blazar jets it can only be a radiatively subdominant component or can dominate only during transient events.

The authors would like to thank the anonymous referee for useful comments and suggestions. The authors also thank A. Mastichiadis, E. Resconi, and M. Huber for useful discussions and comments. I.L. would also like to thank the Department of Astrophysical Sciences at Princeton University for its hospitality during which this project was conceived. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018). M.P. acknowledges support from the Lyman Jr. Spitzer Postdoctoral Fellowship and the Fermi Guest Investigation grant 80NSSC18K1745.