On the minimum jet power of TEV BL Lac objects in the p-$\gamma$ model

TL;DR

This study investigates the minimum jet power required in p-γ models to explain TeV emissions in BL Lac objects, revealing that these powers often surpass black hole Eddington limits, impacting our understanding of their accretion and emission mechanisms.

Contribution

It introduces a method to derive minimum jet powers based on p-γ interaction efficiency and gamma-ray opacity, highlighting that these powers can exceed Eddington luminosities in TeV BL Lac objects.

Findings

Minimum jet powers often exceed Eddington luminosity.

Detection of TeV emission constrains p-γ interaction efficiency.

Implications for black hole accretion modes and emission origins.

Abstract

We study the requirement on the jet power in the conventional p-$\gamma$ models (photopion production and Bethe-Heitler pair production) for TeV BL Lac objects. We select a sample of TeV BL Lac objects whose SEDs are difficult to be explained by the one-zone leptonic model. Based on the relation between the p-$\gamma$ interaction efficiency and the opacity of $\gamma\gamma$ absorption, we find that the detection of TeV emission poses upper limits on the p-$\gamma$ interaction efficiencies in these sources and hence minimum jet powers can be derived accordingly. We find that the obtained minimum jet powers exceed the Eddington luminosity of the supermassive black holes. Implications for the accretion mode of the supermassive black hole in these BL Lac objects and the origin of their TeV emissions are discussed.

Figures

Tables

Table 1. Table 1: Model parameters for SED fitting. We set γ e , max = 10 7 subscript 𝛾 e max superscript 10 7 \gamma_{\rm e,max}=10^{7} for all sources which will not affect our fitting results. The first six objects are HBLs, the last three objects are IBLs.

object	R	B	${\delta }_D$	$q_1$	$q_2$	${\gamma }_{e,\min }$	${\gamma }_{e,\mathrm{break}}$	$L_{e,\mathrm{inj}}$	$E_c$
	(cm)	(Gauss)						(erg/s)	(TeV)
1ES 0229+200	3.98$\times {10}^{15}$	0.43	12.71	1.32	3.95	3.02$\times {10}^2$	5.62$\times {10}^5$	1.55$\times {10}^{41}$	7.33
1ES 0347-121	1.21$\times {10}^{16}$	0.35	28.79	1.83	3.5	3.14$\times {10}^2$	2.82$\times {10}^5$	3.40$\times {10}^{40}$	3.21
1ES 0414+009	7.08$\times {10}^{15}$	0.75	29.13	1.82	3.1	3.11$\times {10}^2$	3.98$\times {10}^4$	5.64$\times {10}^{40}$	1.13
1ES 1101-232	3.16$\times {10}^{15}$	3.60	24.02	1.41	3.35	1.00$\times {10}^1$	4.75$\times {10}^4$	6.89$\times {10}^{40}$	1.11
1ES 1215+303	1.58$\times {10}^{15}$	3.10	15.04	2.24	5.45	3.08$\times {10}^2$	2.51$\times {10}^4$	2.73$\times {10}^{41}$	0.09
1ES 1218+304	5.98$\times {10}^{15}$	0.36	19.81	1.68	3.6	3.00$\times {10}^2$	1.32$\times {10}^5$	1.33$\times {10}^{41}$	1.47
S5 0716+714	2.59$\times {10}^{16}$	0.53	27.97	2.02	3.4	1.00$\times {10}^1$	1.02$\times {10}^4$	1.61$\times {10}^{42}$	0.73
W Comae	2.00$\times {10}^{15}$	2.13	16.36	2.11	3.75	1.00$\times {10}^0$	2.40$\times {10}^4$	2.97$\times {10}^{41}$	0.26
TXS 0506+056	8.91$\times {10}^{16}$	0.11	16.51	1.92	4.07	1.00$\times {10}^1$	1.91$\times {10}^4$	9.21$\times {10}^{42}$	0.38

Table 2. Table 2: Parameters relevant to the jet powers of selected BL Lacs. z 𝑧 z is the redshift of the source; P e , k subscript 𝑃 e k P_{\rm e,k} is the kinetic power in relativistic electrons in the AGN frame in unit of erg/s; Log M BH subscript 𝑀 BH M_{\rm BH} is the logarithm of the SMBH in units of solar masses, M ⊙ subscript 𝑀 direct-product M_{\odot} ; L TeV subscript 𝐿 TeV L_{\rm TeV} is the intrinsic beam-corrected luminosity of the TeV data in units of erg/s; P p , k / L Edd subscript 𝑃 p k subscript 𝐿 Edd P_{\rm p,k}/L_{\rm Edd} is the ratio of the minimum injection proton luminosity to the Eddington luminosity in the case of α 𝛼 \alpha = 10 and q 𝑞 q = 2; P p , k / P e , k subscript 𝑃 p k subscript 𝑃 e k P_{\rm p,k}/P_{\rm e,k} is ratio of the minimum injection proton luminosity to the injection electron luminosity in the case of α 𝛼 \alpha = 10 and q 𝑞 q = 2. For 1ES 0414+009, 1ES 1101-232 and TXS 0506+056, in absence of an estimated black hole mass, we considered an average value of 10 9 ​ M ⊙ superscript 10 9 subscript 𝑀 direct-product 10^{9}M_{\odot} (Paliya et al., 2017 ) .

Equations

$$Q_e({\gamma }_e)=Q_{\mathrm{e},0}{\gamma }_{\mathrm{e}}^{-q_1}{\left[1+{\left(\frac{{\gamma }_e}{{\gamma }_{e,\mathrm{b}}}\right)}^{(q_2-q_1)}\right]}^{-1},{\gamma }_{e,\min }<{\gamma }_e<{\gamma }_{e,\max },$$

$$\begin{array}{r}{N}_{\mathrm{e}}({\gamma }_{\mathrm{e}})=N_{\mathrm{e},0}{\gamma }_{\mathrm{e}}^{-q_1}{\left[1+{\left(\frac{{\gamma }_{\mathrm{e}}}{{\gamma }_{\mathrm{e},\mathrm{b}}}\right)}^{(q_2-q_1)}\right]}^{-1}\\ \times (1+\frac{{\gamma }_{\mathrm{e}}}{{\gamma }_{\mathrm{cool}}}{)}^{-1},{\gamma }_{\mathrm{e},\min }<{\gamma }_{\mathrm{e}}<{\gamma }_{\mathrm{e},\max },\end{array}$$

$$P_{\mathrm{e},\mathrm{k}}\simeq \pi R^2{\Gamma }^2{m}_{\mathrm{e}}{c}^3\int {\gamma }_{\mathrm{e}}{N}_{\mathrm{e}}({\gamma }_{\mathrm{e}})d{\gamma }_{\mathrm{e}}.$$

$$Q_{\mathrm{p}}({\gamma }_{\mathrm{p}})=Q_{\mathrm{p},0}{\gamma }_{\mathrm{p}}^{-q},{\gamma }_{\mathrm{p},\min }<{\gamma }_{\mathrm{p}}<{\gamma }_{\mathrm{p},\max },$$

$$N_{\mathrm{p}}({\gamma }_{\mathrm{p}})\approx t_{\mathrm{dyn}}{Q}_{\mathrm{p}}({\gamma }_{\mathrm{p}}),$$

$$P_{p,k}\simeq \pi R^2{\Gamma }^2{m}_{\mathrm{p}}{c}^3\int {\gamma }_{\mathrm{p}}{N}_{\mathrm{p}}({\gamma }_{\mathrm{p}})d{\gamma }_{\mathrm{p}}=\pi R^3{\Gamma }^2{m}_{\mathrm{p}}{c}^2\int {\gamma }_{\mathrm{p}}{Q}_{\mathrm{p}}({\gamma }_{\mathrm{p}})d{\gamma }_{\mathrm{p}}$$

$$t_{\mathrm{acc}}\simeq \frac{3\alpha }{20}\frac{r_{\mathrm{L}}}{c}\simeq \frac{20\alpha }{3}\frac{{\gamma }_{\mathrm{p}}{m}_{\mathrm{p}}c}{eB}$$

$$t_{\mathrm{esc}}=\frac{R^2}{4D}=\frac{3eB{R}^2}{4\alpha {\gamma }_{\mathrm{p}}{m}_{\mathrm{p}}{c}^3},$$

$${\gamma }_{p,\max }=\sqrt{\frac{80}{9}}\frac{eBR}{\alpha m_p{c}^2}\simeq 10^8{\left(\frac{\alpha }{10}\right)}^{-1}\left(\frac{B}{1\mathrm{G}}\right)\left(\frac{R}{10^{16}\mathrm{cm}}\right)$$

$$f_{\mathrm{ph}}({\gamma }_{\mathrm{p}})=\frac{t_{\mathrm{dyn}}}{t_{\mathrm{ph}}}=\frac{R}{2{\gamma }_{\mathrm{p}}^2}{\int }_{{ϵ}_{\mathrm{th}}^{\mathrm{ph}}/2{\gamma }_{\mathrm{p}}}^{\infty }dϵ\frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^2}{\int }_{{ϵ}_{\mathrm{th}}^{\mathrm{ph}}}^{2ϵ{\gamma }_p}d{ϵ}_r{\sigma }_{\mathrm{ph}}({ϵ}_r)K_{\mathrm{ph}}({ϵ}_r){ϵ}_r$$

$$f_{\mathrm{BH}}({\gamma }_{\mathrm{p}})=\frac{t_{\mathrm{dyn}}}{t_{\mathrm{BH}}}=\frac{R}{2{\gamma }_{\mathrm{p}}^2}{\int }_{{ϵ}_{\mathrm{th}}^{\mathrm{BH}}/2{\gamma }_{\mathrm{p}}}^{\infty }dϵ\frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^2}{\int }_{{ϵ}_{\mathrm{th}}^{\mathrm{BH}}}^{2ϵ{\gamma }_p}d{ϵ}_r{\sigma }_{\mathrm{BH}}({ϵ}_r)K_{\mathrm{BH}}({ϵ}_r){ϵ}_r,$$

$${\tau }_{\gamma \gamma }({ϵ}_1)=\frac{t_{\mathrm{dyn}}}{t_{\gamma \gamma }}=\frac{2R(m_e{c}^2{)}^4}{{ϵ}_1^2}{\int }_{m_e^2{c}^4/{ϵ}_1}^{\infty }dϵ\frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^2}{\int }_1^{\frac{ϵ{ϵ}_1}{m_e^2{c}^4}}ds{\sigma }_{\gamma \gamma }(s)s,$$

$$f_{\mathrm{ph}}(E_{p,c})\simeq 10^{-3}{\tau }_{\gamma \gamma }(E_{\mathrm{c}}).$$

$$L_{\mathrm{ph}}^{\mathrm{UL}}=\frac{4}{3}\pi R^3{\delta }_{\mathrm{D}}^2{m}_{\mathrm{p}}{c}^2\int \frac{5}{8}{f}_{\mathrm{ph}}^{\mathrm{UL}}({\gamma }_{\mathrm{p}}){\gamma }_{\mathrm{p}}{Q}_{\mathrm{p}}({\gamma }_{\mathrm{p}})d{\gamma }_{\mathrm{p}},$$

$$L_{\mathrm{BH}}^{\mathrm{UL}}=\frac{4}{3}\pi R^3{\delta }_{\mathrm{D}}^2{m}_{\mathrm{p}}{c}^2\int f_{\mathrm{BH}}^{\mathrm{UL}}({\gamma }_{\mathrm{p}}){\gamma }_{\mathrm{p}}{Q}_{\mathrm{p}}({\gamma }_{\mathrm{p}})d{\gamma }_{\mathrm{p}}.$$

$$P_{\mathrm{p},\mathrm{k}}\ge L_{\mathrm{TeV}}{\left[\frac{4{\int }_1^{{\gamma }_{p,\max }}\left(\frac{5}{8}{f}_{\mathrm{ph}}+f_{\mathrm{BH}}\right){\gamma }_p^{1-q}d{\gamma }_p}{3{\int }_1^{{\gamma }_{p,\max }}{\gamma }_p^{1-q}d{\gamma }_p}\right]}^{-1}$$

$${\delta }_{\mathrm{D}}>12g^{-2}({\alpha }_1,{\alpha }_2)\left(\frac{B}{1\mathrm{G}}\right){\left(\frac{{\nu }_{\mathrm{syn}}}{10^{18}\mathrm{Hz}}\right)}^{-1}{\left(\frac{E_{\mathrm{c}}}{1\mathrm{TeV}}\right)}^2$$

Full text

galaxies: active — galaxies: jets — radiation mechanisms: non-thermal

1 Introduction

Blazars are the most extreme form of active galactic nuclei (AGN), with their jets pointing in the direction of the observer (Urry & Padovani, 1995). Multi-wavelength observations show that the spectral energy distributions (SEDs) of blazars generally exhibit a two-bump structure. The origin of the low-energy bump is generally considered to be synchrotron radiation of relativistic electrons accelerated in the jet, while the origin of the high-energy bump is still under debate. In leptonic models, the high-energy bump is explained as inverse Compton (IC) scattering, in which the high-energy electrons in the jet up-scatter the low energy photons from the external photon field such as the emission of the broad line region or the accretion disk (external Compton, EC), or the synchrotron radiation of the electrons of the same population (i.e., synchrotron-self Compton, SSC). In hadronic models, the high-energy bump is instead assumed to originate from proton-synchrotron emission, or emission from secondary particles generated in photohadronic and Bethe-Heitler (BH) interactions (Mannheim, 1993; Aharonian, 2000; Atoyan & Dermer, 2003). Hereafter, we denote the photohadronic and BH interactions collectively by $p-\gamma$ interactions.

So far, 71 blazars have been detected in the TeV band, most of which are high-synchrotron-peaked BL Lac objects (HBLs)111http://tevcat.uchicago.edu/. Due to the lack of the strong emission from the external photon field in BL Lac objects, the SSC model is usually employed to explain the high-energy emissions in the leptonic model (Mastichiadis & Kirk, 1997; Petry et al., 2000; Krawczynski et al., 2002; Yan et al., 2014). The TeV emission from blazars is absorbed due to the $\gamma\gamma$ pair production by interacting with the extragalactic background light (EBL). After correcting for EBL absorption, the intrinsic TeV spectrum is harder than the observed one. Particularly, the TeV spectra of some HBLs are too hard to be explained with the SSC mechanism, since the Klein-Nishina (KN) effect softens the IC spectrum in the TeV band. Thus, the hard TeV spectra pose a challenge to the leptonic explanation and may suggest a hadronic origin.

Among hadronic models, the proton synchrotron model is often employed. In general, proton synchrotron spectra typically peak at multi-GeV rather than TeV energies (Böttcher et al., 2013; Paliya et al., 2018), unless one considers extreme scenarios, such as protons with energy $10^{20}\,$eV radiating in kilo-Gauss magnetic fields, or Doppler factors of $\sim 100$ for the jet. Also, some studies (Zdziarski & Böttcher, 2015; Petropoulou & Dermer, 2016) suggest that the minimum jet power in the proton synchrotron model will exceed the Eddington luminosity of the supermassive black hole (SMBH) which launches the jet. In the $p-\gamma$ model, the energy of relativistic protons is mainly lost through photohadronic interactions ($p+\gamma\rightarrow p/n+\pi^{0}+\pi^{\pm}$) and BH pair production ($p+\gamma\rightarrow p+e^{\pm}$), with the radiation zone being relatively compact. However, Sikora et al. (2009) and Sikora (2011) argue that the $p-\gamma$ interactions are very inefficient in flat spectrum radio quasars (FSRQs) so that an extremely high proton power is required in order to explain the high energy radiation with $p-\gamma$ interactions.

It should be noted that there is a robust connection between the efficiency of $p-\gamma$ interactions and the opacity of the internal $\gamma\gamma$ pair production in $p-\gamma$ models, since the target photon fields of these processes are the same. It has been shown that the interaction efficiency of photohadronic processes in the high-energy limit is about 1000 times smaller than the peak $\gamma\gamma$ opacity (Aharonian, 2000; Dermer et al., 2007, 2012). Such a relation implies that if $p-\gamma$ interactions are very efficient, high-energy gamma-ray emission should not be expected to be detected from the same object. On the other hand, the detection of high-energy gamma-ray emission from certain BL Lac objects can in turn place an upper limit on the efficiency of $p-\gamma$ interactions, which then translates to a minimum proton power in the jet. In this work, we will derive conservative yet robust lower limits on jet powers based on observations of some TeV BL Lac objects, utilizing the relation between the $p-\gamma$ interaction efficiency and the internal $\gamma\gamma$ pair production opacity. Note that any emission from $p-\gamma$ interactions, despite the complicated electromagnetic cascade induced by secondary particles, eventually originates from the energy of protons lost in $p-\gamma$ interactions. Thus, the constraint arising from the opacity of the $\gamma\gamma$ annihilation applies to any model in the framework of $p-\gamma$ processes, no matter which radiation mechanism (e.g., synchrotron, IC, pionic radiation) or which type of radiating particles (e.g., electron/positron, muon, neutral pion) are involved.

The rest of this paper is structured as follows. In Section 2 we describe our method to obtain a lower limit on the proton power of an AGN jet. We apply our method to a sample of 9 TeV BL Lacs in Section 3; in Section 4 we present our discussion and conclusions. Throughout the paper, the $\Lambda$CDM cosmology with $H_{\rm{0}}={70\rm{km~{}s^{-1}Mpc^{-1}}}$, $\Omega_{\rm{m}}=0.3$, $\Omega_{\rm{\Lambda}}=0.7$ is adopted.

2 Model description

2.1 The injected particle energy distribution

In BL Lac objects, the target photon field for the $p-\gamma$ interactions and the $\gamma\gamma$ pair production is mainly provided by the synchrotron radiation of electrons. The proton spectral shape is crucial to the overall radiation efficiency of protons. Thus, we firstly model the electron and proton spectrum in the radiation zone of the jet. The parameters are measured in the comoving frame of the jet unless otherwise specified.

We assume a single spherical radiation zone of radius $R$ being composed of a plasma of electrons and protons in a uniformly entangled magnetic field ($B$), and the observed emission is boosted by a relativistic Doppler factor $\delta_{\rm D}$. Assuming the jet moves with a bulk Lorentz factor $\Gamma$ (or with a velocity of $\beta=\sqrt{1-1/\Gamma^{2}}$ in units of the speed of light $c$), we have $\delta_{\rm D}=[\Gamma(1-\beta\rm cos\theta)]^{-1}\approx\Gamma$ for a relativistic jet close to the line of sight in blazars with a viewing angle of $\theta\lesssim 1/\Gamma$. To explain the low-energy bump in the SED, a broken power-law distribution is required for the electron injection spectrum, i.e.

[TABLE]

where $Q_{e,0}$ is the normalization, $\gamma_{e}$ is the electron Lorentz factor, $\gamma_{e,\rm min}$ is the minimum Lorentz factor, $\gamma_{e,\rm max}$ is the maximum electron Lorentz factor, $\gamma_{e,\rm b}$ is the break electron Lorentz factor, $q_{1}$ and $q_{2}$ are the spectral indices below and above $\gamma_{e,\rm b}$. Given an electron injection luminosity $L_{e,\rm inj}$ in the blob, $Q_{e,0}$ can be obtained by $\int Q_{e}\gamma_{e}m_{e}c^{2}d\gamma_{e}=L_{e,\rm inj}/(4/3\pi R^{3})$ where $m_{e}$ is the mass of an electron. We assume a quasi-steady state is reached, and the injection is balanced by radiative cooling and/or particle escape. The number density of the injected electrons in the radiation zone can be obtained by $Q_{e}t_{e}$, where $t_{e}=min\{t_{\rm cool},t_{\rm dyn}\}$. $t_{\rm cool}=\frac{3m_{e}c}{4(U_{\rm B}+\kappa_{\rm KN}U_{\rm ph})\sigma_{\rm T}\gamma_{e}}$ is the radiative cooling time where $U_{\rm B}=\frac{B^{2}}{8\pi}$ is the energy density of the comoving magnetic field, $U_{\rm ph}$ is the energy density of the soft photons, $\sigma_{\rm T}$ is the Thomson scattering cross section and $\kappa_{\rm KN}$ is a numerical factor accounting for the KN effect. $t_{\rm dyn}$ is the dynamical timescale of the blob, which may be determined by the adiabatic expansion of the blob or by the particle injection processes. Typically, we have $t_{\rm dyn}\simeq R/c$. Then, the quasi-steady electron spectrum in the blob shows a double broken power-law form (Inoue & Takahara, 1996), i.e.

[TABLE]

where $N_{\rm e,0}$ is the normalization coefficient which is equal to $Q_{\rm e,0}t_{\rm ad}$, and $\gamma_{\rm cool}=\frac{3m_{e}c^{2}}{4(U_{\rm B}+U_{\rm ph})\sigma_{\rm T}R}$ is the electron Lorentz factor where $t_{\rm ad}=t_{\rm cool}$. Then the kinetic power in relativistic electrons $P_{\rm e,k}$ in the AGN frame is given by

[TABLE]

Protons are assumed to be injected with a power-law distribution222Note that if we assume the proton spectral shape is same as the electron spectral shape, which was a broken power law, the overall $p-\gamma$ interaction efficiency will be decreased, since the amount of protons at high energy, where the $p-\gamma$ interaction efficiency is relatively high, will be reduced significantly. As a result, the required energy budget for relativistic protons will be even higher. Please refer to the discussion in later sections., i.e.

[TABLE]

where $Q_{\rm p,0}$ is the proton injection constant, $\gamma_{\rm p}$ is the proton Lorentz factor, $q$ is the spectral index, $\gamma_{\rm p,min}$ is the minimum proton Lorentz factor which is usually $\sim 1$ and $\gamma_{\rm p,max}$ is the maximum proton Lorentz factor which can be obtained by comparing the acceleration timescale and the escape timescale of protons. For a proton injection rate $Q_{\rm p}(\gamma_{\rm p})$ per unit volume in the blob, the quasi-steady-state proton energy distribution is given by

[TABLE]

since protons generally are not cooled efficiently in the $p-\gamma$ model. Therefore, we can obtain the kinetic power in relativistic protons as

[TABLE]

We assume that particle acceleration is dominated by diffusive shock acceleration, for which the acceleration timescale in the relativistic limit can be evaluated by (Protheroe & Clay, 2004; Rieger et al., 2007)

[TABLE]

under the quasi-linear theory, where $r_{\rm L}$ is the Larmor radius of the proton and $\alpha$ is the ratio of the mean magnetic field energy density to the turbulent magnetic field energy density. Generally, we expect $\alpha=10-100$ or even larger (Lagage & Cesarsky, 1983; Hillas, 1984, but also see e.g. Bell 2004 for the discussion on saturation of turbulent magnetic field), but in the very limiting case the value of $\alpha$ may approach unity. On the other hand, we assume protons escape via diffusion so the escape timescale can be written as

[TABLE]

where $D$ is the diffusion coefficient, $D=\alpha r_{\rm L}c/3$. Of course, there may be other processes which could play a more important role in limiting the acceleration, such as the advective escape or adiabatic cooling. A more sophisticated treatment requires detailed modeling of,e.g., the geometry and the configuration of the magnetic field. We here simply consider diffusion as the main escape mechanism so that we can obtain the maximum proton Lorentz factor as

[TABLE]

2.2 The minimum injection proton power

Due to the large number of free parameters in $p-\gamma$ models, it is generally impossible to obtain a unique set of parameters by modeling the SED. However, it is meaningful to search for the minimum jet power and compare it with the Eddington luminosity of the SMBH. In this subsection, we propose a method of obtaining a robust lower limit on the proton power in the $p-\gamma$ model.

The efficiency of photohadronic interactions and BH (Bethe-Heitler) pair production in a radiation field can be written as (Stecker, 1968; Berezinskii et al., 1990)

[TABLE]

and

[TABLE]

respectively, where $\epsilon$ represents the photon energy in the jet frame, and $\epsilon_{\rm th}^{\rm ph/\rm BH}$ is the photon threshold energy in the rest frame of proton for the photohadronic and BH processes, respectively, $n_{\rm ph}(\epsilon)$ is the number density of the soft photons in the comving frame which is mainly provided by the synchrotron radiation of primary electrons, $\epsilon_{r}$ is the photon energy in the rest frame of proton, $\sigma_{\rm ph}(\epsilon_{r})$ (Mücke et al., 2000) is the cross section for photopion production, $\sigma_{\rm BH}(\epsilon_{r})$ (Chodorowski, 1992) is the cross section for BH pair production, $K_{\rm ph}(\epsilon_{\rm r})$ (Mücke et al., 2000) is the inelasticity of photohadronic interactions and $K_{\rm BH}(\epsilon_{\rm r})$ (Chodorowski, 1992) is the inelasticity of BH pair production.

Generally, the high energy photons produced in the jet will be attenuated by interacting with the synchrotron radiation of primary electrons in the jet. This internal $\gamma\gamma$ absorption optical depth can be calculated as (Finke et al., 2008; Dermer & Menon, 2009)

[TABLE]

where $\sigma_{\rm\gamma\gamma}$ is the $\gamma\gamma$ pair-production cross section, $\sqrt{s}$ is the center-of-momentum frame Lorentz factor of the produced electron and positron (Dermer & Menon, 2009), and $\epsilon_{1}$ is the energy of $\gamma$-ray photons.

To get efficient hadronic emission, we may in principle adjust model parameters to result in a large $f_{p\gamma}(\equiv f_{\rm ph}+f_{\rm BH})$. However, $\tau_{\rm\gamma\gamma}$ will be increased simultaneously since the target photon fields for $\gamma\gamma$ absorption and for $p-\gamma$ processes are the same. Here we define a critical energy $E_{\rm c}$ in the AGN frame beyond which the TeV spectrum of a BL Lac shows a cutoff or softening. If there is no such feature in the TeV spectrum then $E_{\rm c}$ is defined to be equal to the highest energy that the TeV detection extends to. Assuming $\tau_{\rm\gamma\gamma}(E_{\rm c})=1$ for the BL Lac object, we can obtain an upper limit on $f_{p\gamma}$. A simplified expression for this relation can be obtained using the $\delta$-approximation for the cross sections of both processes. In this approximation, the energies of the soft photon $\epsilon_{s}$ and the high energy photon $E_{\rm c}$ in the observer’s frame satisfy the relation $\epsilon_{s}E_{\rm c}\approx 4\delta_{\rm D}^{2}m_{e}^{2}c^{4}$ for the $\gamma\gamma$ annihilation, while the same photons interact with protons of energy $E_{p,c}\simeq\delta_{D}^{2}m_{p}c^{2}\epsilon_{\Delta}/\epsilon_{s}$ where $\epsilon_{\Delta}\simeq 0.3\,$GeV is the energy for the $\Delta$ resonance. Thus, we have $E_{p,c}\simeq 3\times 10^{5}E_{\rm c}$ which does not depend on any model parameters. If we compare the peak cross section for $\gamma\gamma$ annihilation $\sigma_{\rm\gamma\gamma,\rm peak}\simeq 10^{-25}\rm\,cm^{-2}$ with the peak value of the product of the cross section and the inelasticity for the photopion production (i.e.,$\sim 10^{-28}\rm\,cm^{-2}$), we obtain the relation.

[TABLE]

Thus, the condition $\tau_{\rm\gamma\gamma}(E_{\rm c})=1$ suppresses the $p-\gamma$ interaction efficiency to a quite low level, i.e., $f_{\rm ph}(E_{p,c})\sim 10^{-3}$. For protons with other energies, $f_{p\gamma}=f_{p\gamma}(E_{p,c})F(E_{p})$, where $F(E_{p})$ is a normalized function depicting how $f_{p\gamma}$ changes with $E_{p}$. Since we fix $\tau_{\gamma\gamma}(E_{\rm c})=1$ in our treatment, the uncertainty of $f_{p\gamma}$ mainly originates from the uncertainty of $F(E_{p})$ which is determined by the spectral shape of the target photon field for $p-\gamma$ interactions, i.e., the spectral shape of the electron synchrotron radiation. Although the low-energy SED may be fitted with different combinations of parameters, the resulting spectral shape should be always compatible with the observation from the optical band to the X-ray band. Thus, $F(E_{p})$ is more or less fixed for a given source. In addition, the condition $\tau_{\gamma\gamma}(E_{\rm c})=1$ also reduces the degeneracy of the model parameters. Therefore, $F(E_{p})$ is not expected to vary significantly. A similar relation can be also obtained for the BH process. Fig. 1 shows the the $f_{\rm ph}$, $f_{\rm BH}$ and $\tau_{\gamma\gamma}$ as functions of proton/photon energy.

Here we take 1ES 0229+200 as an example to further interpret the relation between $\tau_{\rm\gamma\gamma}$ and $f_{\rm ph}$. The critical energy $E_{\rm c}$ for this source is 7.3 TeV. By adjusting the physical parameters to make $\tau_{\rm\gamma\gamma}(E_{\rm c})$ equal to 1, we have $f_{\rm ph}\simeq 10^{-3}$ at the proton energy $E_{p,c}\simeq 3\times 10^{5}E_{\rm c}\simeq 2.2\times 10^{18}$eV where $p-\gamma$ interaction efficiency reaches the maximum. Since $F(E_{p})$ is determined by the photon spectra of low-energy SED, we can fully determine $f_{p\gamma}=f_{p\gamma}(E_{p,c})F(E_{p})$ for this source, as shown in Fig. 1. One can in principle further increases the value of $f_{p\gamma}$ by adjusting certain parameters, but the gamma-ray opacity around $E_{\rm c}$ will then become larger than unity and we would expect a break or cutoff in the TeV spectra around $E_{\rm c}$, which however, is not seen in the data.

More quantitatively, for a BL Lac object with a hard TeV spectrum, we first correct for the influence of EBL attenuation on the spectrum. Then, we look for the highest energy data point after which a significant suppression or a softening appears in the spectrum, or simply the highest energy data point if no spectrum suppression or softening appears. The corresponding energy of the data point has been defined as $E_{\rm c}$. We then fit the low-energy bump in the SED of the BL Lac object with the synchrotron radiation of primary electrons and choose model parameters to achieve $\tau_{\rm\gamma\gamma}(E_{\rm c})=1$. Based on the resulting model parameters, we calculate the $p-\gamma$ interaction efficiency and denote the obtained value by $f_{p\gamma}^{\rm UL}$, where “UL” means that the obtained $p-\gamma$ efficiency is the upper limit for the source. The maximum beam-corrected luminosity (i.e., assuming that the inferred luminosity is emitted only in a beam of opening angle $\theta_{j}\sim 1/\Gamma$) of electromagnetic particles produced in the $p-\gamma$ interactions can then be given by

[TABLE]

for the photohadronic process where the factor $5/8$ considers about $3/8$ of the lost proton energy goes into neutrinos, and

[TABLE]

for the BH process. We obtain the total beam-corrected luminosity from the $p-\gamma$ processes as $L_{p\gamma}^{\rm UL}=L_{\rm ph}^{\rm UL}+L_{\rm BH}^{\rm UL}$. Regardless of the details of the electromagnetic cascade induced by those electromagnetic particles, we have $L_{\rm TeV}<L_{p\gamma}^{\rm UL}$ simply from the perspective of the energy budget where $L_{\rm TeV}$ is the intrinsic (i.e., beam-corrected) TeV gamma-ray luminosity that cannot be explained by leptonic processes, since the TeV emission eventually originates from the electromagnetic particles generated in the $p-\gamma$ model as we mentioned earlier. On the other hand, according to the expression for the kinetic power in relativistic protons $P_{\rm p,k}$, i.e., Eq. (6), we can see that the ratio of $L_{p\gamma}$ to $P_{\rm p,k}$ does not depend on the proton luminosity, but on the interaction efficiency and the spectral shape of injected protons (i.e., power-law index $q$ and the cutoff energy). Thus, to account for the TeV emission through the $p-\gamma$ processes, the required nonthermal proton injection power of the jet should satisfy

[TABLE]

since $L_{\rm TeV}\leq L_{p\gamma}^{\rm UL}$. We note that $\tau_{\gamma\gamma}(E_{\rm c})$ is not necessarily equal to unity. The cutoff or softening feature in TeV spectrum may be simply due to the cutoff or softening in the spectrum of the emitting particles instead of $\gamma\gamma$ absorption. Therefore, the realistic value of $f_{p\gamma}$ can be even smaller than the one obtained by imposing $\tau_{\gamma\gamma}(E_{\rm c})=1$.

3 Applications

We now apply the procedure introduced in the preceding section to some TeV BL Lac objects. To reduce the uncertainty caused by the model parameters as much as possible, we select our sample of BL Lac objects according to the following criteria: (i) the redshift is known; (ii) (quasi-)simultaneous multi-wavelength SED data are available; (iii) the TeV emission of the sources can not be well reproduced by the leptonic model or at least the origin of the TeV emission is under debate in the previous literature (see the Appendix.A for details). Note that SEDs of these BL Lac objects may still be fitted by the leptonic model (or other models different from $p-\gamma$ models) via introducing geometry effects or multiple emission zones, but we here consider the conventional hadronic scenario and examine the requirement on the jet power in the framework of $p-\gamma$ models.

Based on the above criteria, we collect a sample of 3 IBLs and 6 HBLs (according to the classification by Abdo et al. (2010)) from the TeVCat. The (quasi-)simultaneous SEDs of 1ES 0229+200, 1ES 0347-121, 1ES 0414+009, 1ES 1101-232, 1ES 1215+303, 1ES 1218+304, S5 0716+714, W Comae and TXS 0506+056 are taken from Aliu et al. (2014), Aharonian et al. (2007), H.E.S.S. Collaboration et al. (2012), Aharonian et al. (2007), Aleksić et al. (2012), Rüger, Spanier, & Mannheim (2010), Anderhub et al. (2009), Acciari et al. (2009) and IceCube et al. (2018), respectively. It should be noted that there is no simultaneous GeV data for 1ES 0347-121, 1ES 1101-232, 1ES 1218+304, S5 0716+714 during the TeV observation since the TeV observation is performed before the launch of Fermi satellite. We then use the two-year average Fermi-LAT data (2008.08.04 - 2010.08.01)333The database of the Space Science Data Center (SSDC) SED Builder provides the two-year average Fermi-LAT data; https://tools.ssdc.asi.it/SED/ in our fitting. In addition, the quasi-simultaneous GeV data of W Comae are taken from Abdo et al. (2010). In the SED of 1ES 1218+304, the thermal radiation is prominent in the infrared band which is believed to originate from the host galaxy (Rüger, Spanier, & Mannheim, 2010), and hence we do not consider the interaction processes (including IC radiation, $\gamma\gamma$ annihilation and $p-\gamma$ interactions) on this thermal component. Note that $f_{p\gamma}$ and $\tau_{\gamma\gamma}$ caused by this thermal component still follow the aforementioned relation, so we do not expect that including this external photon field can significantly increase $f_{p\gamma}$ even if this component is from a much more compact region with a larger photon density. We also note that the VHE emission of 1ES 1215+303, 1ES 1218+304, W Comae and TXS 0506+056 show significant variation, while no evidence of variability in TeV is found for the other sources. Our method is based on the relation between the $\gamma\gamma$ annihilation opacity and the $p-\gamma$ interaction efficiency. Thus, the method can be generalized to both high and low state of the blazar as long as these two processes operate in the same target photon field which is considered to be the synchrotron radiation of accelerated electrons.

Using the standard one-zone leptonic SSC model (Katarzyński et al., 2001), we firstly fit the low-energy bump with synchrotron radiation of primary electrons for each BL Lac in the sample. Instead of exploring the entire parameter space to optimize the fitting, we look for parameters to achieve a $\gamma\gamma$ annihilation opacity equal to unity at the critical energy (i.e., $\tau_{\gamma\gamma}(E_{\rm c})=1$ as shown with red data points in Fig. 2), such that the upper limit of the $p-\gamma$ interaction efficiency can be obtained. On the premise of a reasonable fitting to the low-energy bump, we also try to fit the high-energy data with SSC emission as much as possible. The fitting results are shown in Fig. 2 and the model parameters are shown in Table 1. One can see that the SEDs of these 9 BL Lacs from optical to the GeV band can be fitted well by the leptonic model, but it starts to fail for photons above 100 GeV. We then calculate the total luminosity beyond that energy and obtain $L_{\rm TeV}$444The integrated TeV luminosity is estimated by the following methods: Given a series of data points of the energy and the flux $\{(E_{0},F_{0}),(E_{1},F_{1}),\cdot\cdot\cdot,(E_{i},E_{i}),\cdot\cdot\cdot,(x_{n},y_{n})\}$ ($E_{0}<E_{1}<\cdot\cdot\cdot<E_{n}$) which are not fitted by the leptonic model. We calculate the beam-corrected luminosity in the range of $E_{0}-E_{n}$ by the trapezoidal rule, i.e., $L_{\rm TeV}=4\pi D_{L}^{2}\sum_{i=0}^{n-1}(F_{i}+F_{i+1})(E_{i+1}-E_{i})/2\delta_{D}^{2}$ where $D_{L}$ is the luminosity distance of the source. Error bars are not considered..

Considering the obtained leptonic emission as the target photon field, we calculate the $p-\gamma$ interaction efficiency via Eq. 10 and Eq. 11. To study the influence of the proton injection spectrum, we employ several values for the spectral index $q$ in the range of $1.6-2.4$ and three different $\gamma_{\rm{p,max}}$ with $\alpha=1$, $\alpha=10$ and $\alpha=100$ respectively. The resulting minimum $P_{\rm p,k}$ for different BL Lacs in each combination of $q$ and $\gamma_{p,\rm max}$ are shown in Fig. 3. We can see that the minimum $P_{\rm p,k}$ decreases as $q$ and $\alpha$ become smaller. This is because given a harder injection spectrum (i.e., a smaller $q$) and a higher cutoff energy in the spectrum (i.e., a smaller $\alpha$), more energy is distributed to high energy where the $p-\gamma$ interaction efficiency is larger. For most sources, the minimum $P_{\rm p,k}$ is larger than the corresponding Eddington luminosity $L_{\rm Edd}$ for most combinations of $\alpha$ and $q$. Particularly, even with a very hard injection proton spectrum (i.e., $q=1.6$) and with the extreme case of $\alpha=1$, $P_{\rm p,k}$ is still larger than the Eddington luminosity for 1ES 0347-121, 1ES 1218+304, and S5 0716+714, suggesting that a super-Eddington jet luminosity is needed in the $p-\gamma$ model.

We note that the minimum jet powers obtained in this work are conservative. First, a considerable fraction of proton energy lost into EM particles in $p-\gamma$ interactions may be reprocessed into the X-ray or lower energy band via synchrotron radiation of the generated electrons/positrons so that we need a larger $L_{p,\rm inj}$ to account for the TeV emission. Second, we intentionally choose parameters to achieve a gamma-ray opacity equal to unity at the critical energy $E_{\rm c}$ when we fit the low-energy bump in the SED of BL Lacs in the sample. However, as mentioned in the preceding discussion, the opacity can be much smaller than unity and hence the $p-\gamma$ interaction efficiency used in this work is most likely an overestimation. In addition, we do not expect the jet to consist only of relativistic protons. Usually, one would expect the jet to contain more cold protons than relativistic protons and thus the obtained relativistic proton power may only constitute a small part of the jet power.

4 Discussion and conclusions

The Eddington luminosity is obtained by balancing the force of radiation pressure and gravity of an object. Although it is not a strict limit on the luminosity of a black hole, the Eddington luminosity is usually regarded as a reasonable approximation for the maximum jet power of a blazar. Among the radiation models for blazar jets, the leptonic models usually require a sub-Eddington jet power since the radiation efficiency of electrons is high. The low radiation efficiency of protons in the hadronic models (either the $p-\gamma$ model or the proton synchrotron model) obtained in this work and previous studies Sikora (2011); Zdziarski & Böttcher (2015), implies a super-Eddington jet power. Such a jet may be powered by other mechanism such as the Blandford-Znajek mechanism (Blandford & Znajek, 1977) which extracts the spin power of the SMBH or by the super-critical accretion. In the former scenario, however, Zdziarski & Böttcher (2015) pointed out that the magnetic fluxes measured through the radio-core shift effect in some blazars rule out the later mechanism. The latter scenario, i.e., super-critical accretion onto SMBHs, has been studied in various works (e.g. Beloborodov, 1998; Volonteri et al., 2015; Sa̧dowski & Narayan, 2015). Sa̧dowski & Narayan (2015) found that powerful jets with super-Eddington luminosity may be able to launch from the SMBH only under some uncommon conditions (such as in a tidal disruption event). However, even if the condition can be satisfied, the super-critical accretion mode can only last a very small fraction of the lifetime of a SMBH as indicated by Zdziarski & Böttcher (2015), otherwise the growth of the SMBH would be too quick. Thus, such an accretion mode can be only applied to a tiny fraction of blazars. Furthermore, simulation (Sa̧dowski & Narayan, 2016) shows that the radiation in the funnel along the axis is supposed to be super-Eddington (which is the case of BL Lacs) when the accretion is super-critical. However, from the non-detection of the spectral feature of the accretion-disk emission in the SED of the BL Lacs in our sample, we can estimate upper limits of the disk luminosity for these sources to be $10^{44}-10^{45}$erg/s which are sub-Eddington. On the other hand, in the picture of jet/disk symbiosis (Falcke & Biermann, 1995), although it is possible that the accretion power are channelled into the kinetic energy of the jet or the wind than into the disk radiation, the theoretical expectation for the ratio between the jet’s kinetic luminosity and the accretion-disk luminosity is $\lesssim 10$ for a large range of reasonable parameters (Donea & Biermann, 1996; Donea & Protheroe, 2003). It is much smaller than the ratio required in the $p-\gamma$ model for BL Lacs in our sample which is $>10^{3}-10^{4}$. Thus, null detection of the accretion-disk emission from these BL Lacs disfavor super-critical accretion in these sources. Besides, jet powers estimated from radio lobes and X-ray cavities (Merloni & Heinz, 2007; Nemmen et al., 2012) are in conflict with the required super-Eddington jet power at a timescale of $1-10\,$Myr. If the $p-\gamma$ model applies, it probably implies a different picture for the accretion of SMBH in blazars than the one depicted by the standard theory.

Madejski et al. (2016) suggest that the jet power can be reduced significantly by introducing a huge amount of positrons to replace protons (in their case, the number density of positrons is 30 times higher than that of protons) in the jet from the point of view of keeping the neutrality of the jet. We note that although this is a possible solution to some sources, it does not apply to the BL Lacs in our sample. This is because the IC radiation of positrons also suffers the KN suppression and hence cannot explain the hard TeV spectrum (at least in the one-zone model).

To summarize, we obtained a conservative yet robust lower limit on the jet power for TeV BL Lacs for which the standard leptonic model does not work well. The detection of TeV photons from BL Lacs imposed an upper limit for the $\gamma\gamma$ annihilation opacity. Since $p-\gamma$ interactions (including photopion production and the BH process) take place in the same target photon field as $\gamma\gamma$ annihilation, the $p-\gamma$ interaction efficiency is linked with the $\gamma\gamma$ opacity. Based on this relation, we obtained an upper limit for the $p-\gamma$ interaction efficiency which translates to the minimum proton power of the jet if $p-\gamma$ interactions are responsible for the TeV emission from these BL Lacs. By applying this approach to a sample of 9 TeV BL Lacs, we found that the minimum injection proton power is larger than the Eddington luminosity for most combinations of $q$ and $\alpha$. If the Eddington luminosity is the maximum luminosity that a blazar jet can achieve, the $p-\gamma$ process may not be responsible for the TeV emission in these TeV BL Lacs. One then may have to consider the leptonic origin with more complicated topology of the radiation zone. On the other hand, the radiation efficiency of protons in the hadronic model that employs the hadronuclear interaction is not related to the gamma-ray opacity (e.g. Dar & Laor, 1997; Bednarek & Protheroe, 1997; Araudo et al., 2010; Barkov et al., 2010; Khangulyan et al., 2013; Liu et al., 2018; Sahakyan, 2018), and hence may provide a solution to fit the TeV spectrum with a sub-Eddington jet power.

Acknowledgments.

Appendix A Difficulties of the standard leptonic models in fitting TeV BL Lacs

We here summarize the difficulties of the leptonic interpretation of the TeV emissions from the sources in our selected sample, based on results and arguments in previous literature:

(1) A large value of the Doppler factor ($\delta_{\rm D}$) is needed. For 1ES 0229+200, one of the important outcome of the SSC interpretation (Aliu et al., 2014) is that the minimum $\delta_{\rm D}$ required in the fitting of SSC model is $\delta_{\rm D}\geqslant 53$, which is significantly higher than the commonly adopted value for blazar jets. In earlier studies, a large $\delta_{\rm D}$ has been commonly suggested for this object, e.g., Tavecchio et al. (2009) adopt $\delta_{\rm D}$ = 50 and Kaufmann et al. (2011) adopt $\delta_{\rm D}$ = 40. Similarly, for 1ES 0414+009, H.E.S.S. Collaboration et al. (2012) adopt $\delta_{\rm D}=40$ with either a SSC model and an EC model. Even with such a large $\delta_{\rm D}$, the SED fitting is only marginal. For 1ES 1215+303, $\delta_{\rm D}=60$ is required (Aleksić et al., 2012). For 1ES 1218+304, an extreme value of $\delta_{\rm D}$ = 80 is employed (Rüger, Spanier, & Mannheim, 2010). An extremely high $\delta_{\rm D}$ imply a very fast movement of the radiation region, though such a high $\delta_{\rm D}$ is inconsistent with radio observations of the movement of knots (Lister et al., 2013) or with the statistics of blazars (Henri & Saugé, 2006).

(2) A high value of $\gamma_{\rm e,min}$ is needed ($\sim 10^{4}$) in the SED fittings for 1ES 0229+200, 1ES 0347-121, 0414+009, 1ES 1101-232 and 1ES 1215+303. It requires some specific conditions in those sources to make such fine-tuned values of $\gamma_{\rm e,min}$ physically reasonable (Katarzyński et al., 2006).

(3) The Fermi-LAT observation is hard to reconcile with the VHE observation under the leptonic model for some sources. For S5 0716+714, either fitting with the one-zone leptonic model or with the spine-sheath model, the predicted gamma-ray flux overshoot the observed flux by Fermi-LAT. Moreover, an extremely high $\gamma_{\rm e,min}=10^{4}$ is also needed (Anderhub et al., 2009). For W Comae, SED was initially well fitted with both the SSC model and the EC model by Acciari et al. (2009) without the Fermi observation. However, the leptonic model fails to fit SED after the GeV data (Böttcher et al., 2013) is included.

(4) If a high-energy neutrino event is coincident both temporally and spatially with a $\gamma$-ray flare from a blazar such as in the case of TXS 0506+056 (IceCube Collaboration et al., 2018), hadronic processes has to be considered. Therefore, TXS 0506+056 is included in our sample.

We here take 1ES 0229+200 as an example to show why the leptonic model can not fit the hard TeV spectrum. The synchrotron peak frequency of electrons can be estimated as $\nu_{\rm syn}\approx 3\times 10^{6}\gamma_{\rm e,b}^{2}B\delta_{\rm D}/(1+z)$. Meanwhile, the SSC spectrum peaks at $\nu_{\rm ssc}\approx 4\gamma_{\rm e,b}^{2}\nu_{\rm syn}/3$ which is emitted by electrons of $\gamma_{\rm e,b}$ scattering off the photons from the peak of the synchrotron bump, if we do not consider the KN effect. Thus, in principle one can ascribe the hard TeV spectrum to the SSC emission as long as $h\nu_{\rm ssc}\geq E_{\rm c}$ which is generally true for IBL and HBL. But one also need to guarantee that KN effect does not interfere the spectrum below $E_{\rm c}$, i.e., $\nu_{\rm KN}\simeq\delta_{\rm D}\gamma_{\rm e,b}g(\alpha_{1},\alpha_{2})m_{\rm e}c^{2}/h(1+z)>E_{\rm c}$ where $g(\alpha_{1},\alpha_{2})={\rm exp}\left[1/(\alpha_{1}-1)+1/2(\alpha_{2}-\alpha_{1}))\right]\lesssim 1$ (Tavecchio et al., 1998) with $\alpha_{1}$ and $\alpha_{2}$ the spectrum index (i.e., $f_{\nu}\propto\nu^{-\alpha_{1,2}}$) of the synchrotron emission below and above the peak respectively. It translates to

[TABLE]

Take 1ES 0229+200 for example, the observational synchrotron peak frequency is $\nu_{\rm syn}\approx 10^{18}$Hz, $E_{\rm c}\simeq 7\,$TeV and $g(\alpha_{1},\alpha_{2})\simeq 0.3$. Thus, to fulfill the above relation, we need $\delta_{\rm D}>5000(B/1\rm\,G)$. If we consider the typical magnetic field $B=0.1-1\,$G, the required $\delta_{\rm D}$ is $>100-1000$ which far exceeds the typical value. If we want to fit the TeV flux with the SSC radiation, we need to impose $\nu_{\rm KN}$ to the TeV range. A large $\gamma_{\rm e,b}$ has to be assumed, however, at the expanse of a poor fitting to the low-energy bump (as shown in Fig. 4).