A Two-Zone Model as origin of Hard TeV Spectrum in Extreme BL Lacs

Near the SMBH, MeV-photons annihilate with lower energy photons ( i.e., $\gamma+\gamma\rightarrow e^{\pm}$), producing an $e^{\pm}$ outflow that moves with velocity $\beta_{\rm pl}\approx 0.3-0.7$. These pairs produce an optically thick environment if the luminosity above $\rm 511\,keV$ is greater than $L_{\rm keV}\gtrsim 3\times 10^{-3}\,L_{\rm Edd}$, where $L_{\rm Edd}\approx 1.26\times 10^{47}\,{\rm erg\,s^{-1}}\,(M_{\bullet}/10^{9}M_{\odot})$ is the Eddington luminosity and $M_{\bullet}$ is the mass of the SMBH. Radiation only escapes at the photosphere’s radius, $R_{\rm ph}$, producing a line shape peaking at $\varepsilon_{\rm pl}\approx 511\,\rm keV$ (Beloborodov, 1999). The photosphere occurs at a compact radius of the order of Schwarzschild radius $R_{\rm g}\approx 1.5\times 10^{14}\,{\rm cm}\,(M_{\bullet}/10^{9}M_{\odot})$, and the released radiation is very anisotropic where the greater contribution is confined within a solid angle $\Omega_{\rm pl}\lesssim 0.2\pi$. Here, we assume the emission’s shape is very narrow such that a delta approximation can describe its photon distribution as (Fraija et al., 2020)

where the normalization constant is $u_{\rm pl}=L_{\rm keV}/(\Omega_{\rm pl}R_{\rm ph}^{2}\beta_{\rm pl}c)$ with $c$ de speed of light.

Photons coming from the BLR have contributions mainly from two emission lines corresponding to H I $\rm Ly\alpha$ and He I $\rm Ly\alpha$ with energies peaking at $E_{\rm BLR}=10.2\,{\rm eV}$ and $40.2\,{\rm eV}$, respectively. This emission, spatially distributed in a thin shell geometry, is thought to the result of the scattering of UV photons from the accretion disk into clouds located at a distance of (e.g. Ghisellini & Tavecchio, 2009; Kaspi et al., 2007)

For BL Lacs, the bolometric BLR’s luminosity has an upper limit given by the relation $L_{\rm BLR}\lesssim 5\times 10^{-4}L_{\rm Edd}$ (Ghisellini et al., 2011), and if this is a fraction of the disc luminosity $L_{\rm BLR}=\phi_{\rm BL}L_{d}$, then the BLR energy density becomes (Hayashida et al., 2012)

The photon distribution of the emission’s line in the delta-approximation is given by

Photons coming from the DT are thought of as reprocessing radiation from the accretion disc. Therefore, we can calculate its luminosity as $L_{\rm DT}=\phi_{\rm DT}L_{d}$ (e.g., Ghisellini & Tavecchio, 2008). The radius of this region can be calculated using the result of Cleary et al. (2007), which is given by

where the emission can be modeled as a grey-body with a temperature of 100-1000 Kelvin and a frequency peak at $\nu_{\rm DT}=2\times 10^{13}\,{\rm Hz}$. If the above relation holds for any blazar, then we observe that the photon energy density is independent of the disc luminosity (Hayashida et al., 2012)

where the photon distribution is normalized by $u_{\rm DT}=\int dE\,E\,n_{\rm DT}(E)$.

In a magnetic field, charged particles move along the field lines losing energy via synchrotron radiation. The power radiated by an electron distribution is (Blumenthal & Gould, 1970; Finke et al., 2008; Sauge & Henri, 2004)

where $e$ is the electron charge, $\hbar$ is the Plank’s constant, $x=\varepsilon^{\prime\prime}/\varepsilon^{\prime\prime}_{c}$, the characteristic energy is given by $\varepsilon^{\prime\prime}_{c}=\frac{3eB^{\prime\prime}\hbar}{2m_{e}c}{\gamma^{\prime\prime}_{e}}^{2}\sin{\theta}$ and ${\theta}$ is the pitch angle, the function $R_{s}$ is given in Finke et al. (2008).

An electron moving inside a radiation field produces Compton scattering. As well, synchrotron electron energy loss for each collision produces high-energy photons. The total emissivity coefficient produced by an isotropic electron population is (Blumenthal & Gould, 1970)

Finally, the IC radiation is calculated using the above equation but with a different photon spectrum for SSC and EC emission. The corresponding synchrotron emission for SSC is $n^{\prime\prime}_{s}(\varepsilon^{\prime\prime}_{s})\approx J^{\prime\prime}_{s}(\varepsilon^{\prime\prime})/(4\pi{R_{o}^{\prime\prime}}^{2}c\varepsilon^{\prime\prime}_{s})$ and for EC is the DT radiation with a greybody type distribution.

The presence of a radiation field produces interactions with accelerated protons via photopion process. The main contribution comes from $\Delta$-resonance, $p+\gamma\rightarrow p+\pi^{0}(n+\pi^{+})$. The final stable particles resulting from pion decay are $\pi^{0}\rightarrow 2\gamma$ and $\pi^{+}\rightarrow\nu_{\mu}+\bar{\nu}_{\mu}+\nu_{e}+e^{+}$. From relativistic kinematics, the energy threshold to pion production in the head-on collision with pair-plasma’s photons is

where $m_{\pi}$ is the pion mass, $m_{p}$ is proton mass and $\varepsilon^{\prime}$ is the photon target in the comoving frame.

The production rate of stable particles is given by (Kelner & Aharonian, 2008)

with $\eta=\frac{4\varepsilon^{\prime}_{p}\varepsilon^{\prime}}{m_{p}^{2}c^{4}}$ and $x=\frac{\varepsilon^{\prime}_{i}}{\varepsilon^{\prime}_{p}}$. The function $\Phi$ contains the physic of the process and depends on the type of particle to be considered. The label ${\rm i}$ represents $\gamma,\nu,\bar{\nu},e^{-},e^{+}$. The energy loss timescale is calculated considering all the particle type interactions as (Kelner & Aharonian, 2008)

with $\varepsilon^{\prime}_{min}=\frac{\eta_{0}m_{p}^{2}c^{4}}{4\varepsilon^{\prime}_{p}}$ and $\eta_{0}\approx 0.303$.

Another mechanism that the proton interacts with the radiation field is via photopair production, $p+\gamma\rightarrow p+e^{-}+e^{+}$ , i.e, Bethe-Heitler process. From relativistic kinematics, The energy threshold to pion production in the head-on collision is

and the electron production rate for photopair process in the case of $\delta$-approximation, assuming an isotropic photon distribution is given by (Romero & Vila, 2008)

where the mean energy’s fraction carried by an electron from the parent proton is approximated as $x\approx m_{e}/m_{p}$ (i.e, $\gamma^{\prime}_{e}\approx\gamma^{\prime}_{p}$) , $\alpha_{f}$ is the fine structure constant, $\kappa=\frac{2\gamma^{\prime}_{p}\varepsilon^{\prime}}{m_{e}}$. Finally, the energy loss timescale is

The functions $\psi$ and $\varphi$ are parametrized following Chodorowski et al. (1992).

Additionally to photohadronic processes, protons radiate via synchrotron if the magnetic field strength is enough to have a significant flux contribution. The flux from proton synchrotron is calculated using the equation (10) but replacing $m_{e}$ by $m_{p}$ (Dermer & Menon, 2009).

Annihilation-line photons will attenuate the gamma-ray spectrum produced in the inner region from the pair-plasma and the radiation produced in the outer blob, BLR, and DT when these are passing through them. The optical depth for $\gamma\gamma$ attenuation is calculated following (Stecker et al., 1992; Dermer & Menon, 2009)

where $\sigma_{\gamma\gamma}$ is the cross-section, $\mu$ the collision angle, and $l$ is the mean distance traveled by photons. Depending on the seed-photon source, the term $l$ corresponds to the size of the inner blob, the width of the shell of the BLR, and the size of the DT. The threshold energy for pair creations and the created pair’s velocity in the center of the mass frame is given by

respectively.

We set the photosphere’s velocity as $\beta_{\rm pl,ph}\approx 0.3$ (Beloborodov, 1999) and the emission line peaking at $\varepsilon_{\rm pl}\approx 511\,{\rm keV}$ with a photon distribution given by equation (3). We expect the formation of an outflowing pair-plasma when the luminosity is $L_{\rm keV}\gtrsim 3\times 10^{-3}L_{\rm Edd}$, as previously mentioned. In addition, as the pair-plasma will be produced by the radiation of the accretion disc, we assume $L_{\rm keV}\sim L_{d}$. However, if the ratio of BLR and disc luminosities, as generally accepted, is $L_{d}=10\,L_{\rm BLR}$, then the disc luminosity cannot be greater than $L_{\rm d}\lesssim 5\times 10^{-3}\,L_{\rm Edd}$ following the classification suggested by Ghisellini et al. (2011). Therefore, we take $L_{d}=5\times 10^{-3}L_{\rm Edd}$ as benchmark value for our sample. This corresponds to a photon energy density of $u_{\rm pl}\simeq 3\times 10^{6}\,(M_{\bullet}/10^{9}M_{\odot})\,\rm erg\,cm^{-3}$ in the pair-plasma frame.

Since this region is expected to be located close to the SMBH, we consider a compact size of a few $R_{g}$. Note that this scale can be associated with a minimum observed variability imposed by causality condition $t_{\rm var}^{\rm ob}\gtrsim r_{g}/\mathcal{D}_{i}\sim\rm few\;minutes$.²²2Rapid variability timescales of this order have been observed in some HBL, e.g., Aharonian et al. (2007); Albert et al. (2007). The variability timescale in the VHE gamma-ray band is not established for the case of TBLs, except for 1ES 0229 +200 and 1ES 1218 +304; the former shows an indication of timescales of years and the last one on days. In that case, our model can explain the observation without violating the causality constraint even in the case of compact size of the order of $R^{\prime}\sim 10^{14}\,\rm cm\lesssim\mathcal{D}_{i}t_{\rm daily}$.

The hadronic contribution from this blob is calculated assuming a proton distribution with a hard spectral index of $\alpha_{p}=1.8$, being this value in the range of predicted by shock acceleration (Stecker et al., 2007). Also, we assume the minimum and maximum proton energies in the comoving frame as $\varepsilon^{\prime}_{p,\rm min}=1\,\rm GeV$ and $\varepsilon^{\prime}_{p,\rm max}=100\,\rm TeV$, respectively. Here, we do not assume a larger value of proton energy because the photons coming from the pair-plasma only interact via photopion process with protons up to that value, although we use a Lorentz factor of $\Gamma_{i}=10$ (see Figure 2a). Also, this Figure shows that secondary electrons produced by the photopair process can be neglected compared with the photopion process contribution because the energy loss timescale is approximately two orders of magnitude greater. As observed in Figure 2a), the resulting spectrum will be a narrow shape that shifts the peak to higher energies and decreases the timescale depending on $\Gamma_{i}$. As mentioned, we chose a hard proton distribution, then $\Gamma_{i}$ cannot be higher than 3 because higher values imply a $\gamma$-rays peak beyond 10 TeV. This statement is easy to see from the approximation $\varepsilon_{\gamma}=0.1\varepsilon_{p}$. As the efficiency is approximately constant in a specific proton energy range, the peak will be near that range’s maximum proton energy value. For example, for a high value of $\Gamma_{i}$ we have $E^{\rm ob}_{\gamma}\simeq 100\,{\rm TeV}\,(\varepsilon^{\prime}_{p}/100\,{\rm TeV})(\Gamma_{i}/10)$, and for a low value we have $E^{\rm ob}_{\gamma}\simeq 10\,{\rm TeV}\,(\varepsilon^{\prime}_{p}/30\,{\rm TeV})(\Gamma_{i}/3)$. From Table 1, we note that 1ES 0229+200 is the only TBL whose expected peak is higher than 10 TeV, whereas the rest could be at lower energies, implying that $\Gamma_{i}$ could be less than 3. Therefore, we expect the Lorentz factor lies in the range of $1.5\leq\Gamma_{i}\leq 3$ for TBL sources. It is worth noting that the photopion process involves neutrino emission via charged pion decay products. Although this topic is beyond the scope of the current manuscript, a simple estimation could be computed by using the gamma-ray luminosity ($L_{\nu}\simeq\frac{1}{3}L_{\gamma}$) and spectral break ($\varepsilon_{\nu}\simeq 0.5\varepsilon_{\gamma}$) ( e.g., Murase et al., 2014; Fraija et al., 2018). For example, for the case of 1ES 0229 +200 and using the values of table 1, we obtain $L_{\nu}\approx 5\times 10^{44}\,\rm erg\,s^{-1}$ and $\varepsilon_{\nu}\approx 6\,\rm TeV$. In a forthcoming paper, we will present a detailed study of neutrino emission of all TBLs. Here, we only present the neutrino flux of 1ES 0229 +200 as representative case (see the brown line of the left-hand top panel in Figure 3). Furthermore, gamma-rays emitted in this region will be attenuated by the radiation of the annihilation line, BLR, outer-blob, and DT before being observed. Figure 2b shows that the whole emission in the range of 1 MeV-10 GeV for $1.5\leq\Gamma_{i}\leq 3$ is strongly attenuated by the annihilation line emission; therefore, the cascade emission of secondary pairs would be suppressed. Additionally, it is essential to mention that although the BLR emission can be reflected in the inner blob spectrum, it has little attenuation effect, mainly at 100 GeV for $M_{\bullet}=10^{9}M_{\odot}$. This effect would be smaller by a factor $\sim 3$ for $M_{\bullet}=10^{8}M_{\odot}$ compared to one with $M_{\bullet}=10^{9}M_{\odot}$ due to the scales of the optical depth $\tau_{\gamma\gamma}\propto r_{\rm BLR}$ and the radius $r_{\rm BLR}\propto L_{d}^{0.5}\propto L_{\rm Edd}^{0.5}\propto M_{\bullet}^{0.5}$. Therefore, only in the case of 1ES 0229+200 the effect is noticeable below TeV energies and then gamma-rays above $1\,\rm TeV$ emerge without significant attenuation. It is worth noting that the DT does not play an attenuation effect at the highest energies.

The strength of the magnetic field in the inner blob could be estimated assuming that the magnetic luminosity is conserved along the jet, which is supported by observations (e.g., see O’Sullivan & Gabuzda, 2009). Furthermore, Konigl (1981) suggested that a magnetic field varying as $B\propto r^{-1}$ can explain the electromagnetic emission of most compact VLBI jets. Here, we use the following expression

where the parameter values of the outer blob will be discussed in next sections. Table 2 shows that our model demands values of $B_{\rm o}\approx(0.1-0.3)\rm\,G$ and $r_{\rm o}\approx 10^{17}\rm cm$, for a daily variability timescale. Then using these values and a distance of the inner blob of $r_{i}\approx 10^{14}\rm\,cm$, the magnetic field in the inner blob becomes $B_{i}\sim 10^{2}\rm\,G$. Accelerated protons confined by an intense magnetic field in this blob cool down by synchrotron process. Protons accelerated at ultra-high energies (a few $\sim{\rm EeV}$) might radiate photons with energies at a few of $\sim{\rm GeV}$, which are strongly attenuated by the radiation of the pair-plasma (see Figure 2). We do not consider the treatment of the proton synchrotron scenario because its effect is small compared with the radiation of primary electrons. On the other hand, the electron distribution in the inner blob is assumed to be described by a power-law function with the same spectral index as the proton one, $\alpha_{e}=\alpha_{p}$ (Abdo et al., 2011). Therefore, assuming the charge neutrality condition, the electron normalization constant $K_{e}$ can be estimated once we determine the proton content. Electrons are cooling mainly by adiabatic expansion and synchrotron losses, and only the second mechanism produces a steeper distribution such that $\alpha_{e,2}=\alpha_{e,1}+1\approx 3$. The condition of equality between adiabatic ($t_{\rm ad}=2R^{\prime}/c$) and synchrotron losses, ($t_{\rm syn}=6\pi m_{e}c/(\sigma_{\rm T}\,B_{i}^{\prime 2}\,\gamma_{e}^{\prime}$) could estimate the energy break of the electron population as (Cerruti et al., 2015)

Note that there is a strong dependency with the value of $B_{i}$, and if this value increases, then the energy break of the electron population shifts to lower energies and viceversa. For instance, if the magnetic field has a value of hundreds of Gauss and electrons are confined inside small regions, the energy break must appear only at relativistic energies avoiding the ultra-relativistic regime $\gamma_{e,\rm br}^{\prime}\gg 10$. Therefore, a strong magnetic field implies a steeper electron distribution, suppressing the synchrotron radiation. The observed energy peak of the electron synchrotron in this blob is determined by the value of the energy break, which is given by the above relation. Then, we obtain the expression

We note that the synchrotron energy peak is very dependent on the strength of the magnetic field, with $B_{i}\sim 100\rm\,G$, and the peak appears in the radio band. Therefore, the data reported for some TBLs in the radio band can help us to constrain the value of the magnetic field via synchrotron radiation in this region

For example, considering the case of 1ES 0229+200 with $\alpha_{e,1}=2$ and $K_{e}^{\prime}$ (determined by the neutrality condition), we obtain

which agrees with the suggested one by the equation (22).

We assume this blob is moving faster than the inner blob with $\Gamma_{\rm o}=5$, and the blob’s size is associated with a variability timescale of one day, corresponding to $R_{\rm o}^{\prime\prime}\lesssim 2.6\times 10^{16}\,{\rm cm}\,(t_{\rm var}/{\rm day})(\mathcal{D}_{\rm o}/10)$. If the blob’s radius corresponds to the jet cross-section and the jet half-opening angle is $\theta_{j}\simeq 1/(\Gamma_{\rm o})$, then the dissipation distance could be estimated by the expression

Moreover, using the equations (4) and (7) the location of the BLR and the DT is $r_{\rm BLR}\approx 7(2)\times 10^{16}{\rm cm}$ and $r_{\rm DT}\approx 1.7(0.6)\times 10^{18}\,{\rm cm}$ for $M_{\bullet}=10^{9}(10^{8})\,M_{\odot}$ and $L_{d}=5\times 10^{-3}L_{\rm Edd}$, respectively. Note that the dissipation region is located outside the BLR, but still immersed in the DT radiation field $r_{\rm BLR}\leq r_{\rm o}\leq r_{\rm DT}$. Therefore, as the outer blob is beyond the BLR influence, its radiation can be discarded, and only the infrared DT photons will be considered in our calculation.