On the Direct Correlation between Gamma-Rays and PeV Neutrinos from Blazars

We use a one-zone model consisting of an isotropic, homogeneous, and spherical emission region, or blob, with radius ${R}_{\mathrm{blob}}^{\prime }$, that moves relativistically with Doppler factor ${{\rm{\Gamma }}}_{\mathrm{bulk}}$. Electrons and protons are injected with power-law spectra, ${d}^{2}n^{\prime} /d\gamma ^{\prime} {dt}^{\prime} =K^{\prime} {\gamma }^{\prime \alpha }$ for ${\gamma }_{\min }^{\prime }\lt \gamma ^{\prime} \lt {\gamma }_{\max }^{\prime }$, where ${d}^{2}n^{\prime} /d\gamma ^{\prime} {dt}^{\prime}$ is the differential particle injection rate per volume, α the power-law index, ${\gamma }_{\min }^{\prime }$ and ${\gamma }_{\max }^{\prime }$ the minimum and maximum Lorentz factor of the particles, and $K^{\prime}$ a normalization factor determined by the particle injection luminosity, ${L}_{\mathrm{inj}}^{\prime }$. The injected electrons are henceforth referred to as primary electrons, whereas we denote as secondary electrons those created by hadronic interactions or $\gamma \gamma$ pair production. We allow primary electrons and protons to have separate parameter values for ${L}_{\mathrm{inj}}^{\prime }$, ${\gamma }_{\min }^{\prime }$, ${\gamma }_{\max }^{\prime }$, and α. The emission region is assumed to be filled with a homogeneous, randomly oriented magnetic field of strength $B^{\prime}$. Neutrinos and "optically thin" photons can freely stream out of the blob on the timescale ${t}_{\mathrm{fs}}^{\prime }=3R^{\prime} /4c$. For simplicity, we assume the escape rate for charged particles, in the slow cooling case, to be a fixed multiple of the free-streaming timescale, ${t}_{\mathrm{esc}}={t}_{\mathrm{fs}}/{f}_{\mathrm{esc}}\,=10\,{t}_{\mathrm{fs}}$.³

The list of relevant interactions includes synchrotron emission and synchrotron self-absorption, IC scattering by both electrons and protons, $\gamma \gamma$ pair production and annihilation, Bethe–Heitler (BH) photo-pair production $p+\gamma \to p\,+{e}^{\pm }$, and photohadronic ($p\gamma$) interactions ($X+\gamma \to X^{\prime} +\pi$), where X and $X^{\prime}$ denote either a proton or a neutron and π includes charged or neutral pions. Secondary particles such as ${\pi }^{\pm }$ and ${\mu }^{\pm }$ can in principle radiate before they decay, but for the parameter values relevant to this study, the effect is negligible. The details are listed in Table 4 and described in Appendix B.

The low-frequency hump in the SED of a generic blazar extends from the radio band to the UV, and in some cases even to the X-ray band; the high-frequency component can be observed from X-rays up to TeV γ-rays. Here we discuss four scenarios that may in principle account for the shape of the SED: the pure leptonic (SSC) model, the lepto-hadronic synchrotron-self-Compton (LH-SSC) model, the lepto-hadronic pion (LHπ) model, the lepto-hadronic proton-synchrotron model, and the proton model, which is a purely hadronic model. The defining features of these models are summarized in Table 1.

Table 1.  List of Models

	First Peak	Middle Range	Second Peak
	(eV–keV)	(keV–MeV)	(MeV–TeV)
SSC	L	L	L
(Pure leptonic)	Primary ${e}^{-}$ synchrotron	SSC	SSC
LH-SSC	L	H	L
(Lepto-hadronic)	Primary ${e}^{-}$ synchrotron	Secondary leptonic	SSC by primary ${e}^{-}$
LH-${\boldsymbol{\pi }}$	L	H	H
(Lepto-hadronic)	Primary ${e}^{-}$ synchrotron	Secondary leptonic	Secondary leptonic or γ-rays from direct ${\pi }^{0}$ decay
LH-psyn	L	H	H
(Proton synchrotron)	Primary ${e}^{-}$ synchrotron	Proton synchrotron or secondary leptonic	Proton synchrotron
Proton	H	H	H
(Pure hadronic)	Proton synchrotron	Secondary leptonic	Secondary leptonic or γ-rays from direct ${\pi }^{0}$ decay

Note. In the table, "L"—"Leptonic" and "H"—"Hadronic." "LH" is the abbreviation for "Lepto-hadronic," which is a mixture of Leptonic and Hadronic components.

Download table as:  ASCIITypeset image

In both the SSC and LH-SSC models, the first hump is described by synchrotron emission of primary electrons, and the high-frequency component is due to IC scattering of those photons by the same electrons. The LH-SSC model contains an additional hadronic component compared to the pure leptonic SSC model, which may fill the gap between the two humps (e.g., accounting for the X-ray emission from PKS B1424-418).

In the "lepto-hadronic" models internally, the low-frequency hump is likewise described by synchrotron emission of primary electrons, whereas the second hump arises from hadronic processes. Depending on the parameters of the source, the dominant contribution to the second peak can be γ-rays from ${\pi }^{0}$ decays, synchrotron and IC radiation emitted by ${e}^{\pm }$ from ${\pi }^{\pm }$ decays, internal $\gamma \gamma$ annihilation, or proton-pair production (BH), where it can be called an "$\mathrm{LH} \mbox{-} \pi$" model. The second peak can also be produced in some cases by the proton-synchrotron model (e.g., FSRQ 3C 279; Diltz et al. 2015), namely, the "LH-psyn" model.

In the proton model, leptonic emission is subdominant at all wavebands. The first peak in the SED is attributed to proton-synchrotron radiation, and the second peak is produced by the same type of hadronic processes in the LH-π model.

While both the leptonic and the hadronic models have successfully explained the SED of a number of blazars, for example Mrk 421, 3C 279, etc., the unique combination of the SED with the PeV neutrino information places stringent limits on the model of PKS B1424-418. Here we use analytical and semi-analytical calculations to demonstrate that neither the lepto-hadronic (LH-π, LH-psyn) nor the purely hadronic proton model can simultaneously explain the SED and the PeV neutrino event, leaving the SSC and LH-SSC models of the SED as the only viable contenders. Analytical arguments also give useful constraints on the parameter space.

Proton model: In pγ interactions, the neutrino energy is roughly 5% that of the parent proton, ${E}_{\nu }\sim 0.05\,{E}_{p}$. For PKS B1424-418 at the redshift z = 1.522, the Lorentz factor of the proton in the comoving frame can be written in terms of the bulk Lorentz factor of the blob Γ and the neutrino energy in the observer frame ${E}_{\nu ,\mathrm{PeV}}^{\mathrm{ob}}\equiv {E}_{\nu }^{\mathrm{ob}}/\mathrm{PeV}$ as

$$\begin{eqnarray}&&{\gamma }_{p}^{\prime }\sim 5\times {10}^{7}\,{{\rm{\Gamma }}}^{-1}\,{{\rm{E}}}_{\nu ,\mathrm{PeV}}^{\mathrm{ob}}.\end{eqnarray} \tag{ 1 }$$

If the low-energy peak is attributed to proton-synchrotron emission, the peak frequency of the synchrotron emission, ${\nu }_{\mathrm{pk},1}$, obeys

$$\begin{eqnarray}&&h{\nu }_{\mathrm{pk},1}^{\prime }={m}_{e}{c}^{2}\displaystyle \frac{{m}_{e}}{{m}_{p}}\left(\displaystyle \frac{B^{\prime} }{{B}_{\mathrm{crit}}}\right){\gamma }_{p}{^{\prime} }^{2},\end{eqnarray} \tag{ 2 }$$

where ${B}_{\mathrm{crit}}\equiv 4.41\times {10}^{13}$ G is the critical magnetic field. Combining Equations (1) and (2) provides a constraint on the magnetic field,

$$\begin{eqnarray}&&B^{\prime} {| }_{\mathrm{pk},1}=(7\times {10}^{-4}\ {\rm{G}})\left(\displaystyle \frac{{\nu }_{\mathrm{pk},1}^{\mathrm{ob}}}{{10}^{14}\ \mathrm{Hz}}\right)\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{bulk}}}{10}\right){\left(\displaystyle \frac{{E}_{\nu }^{\mathrm{ob}}}{\mathrm{PeV}}\right)}^{-2}.\end{eqnarray} \tag{ 3 }$$

The peak frequency of the second hump in the SED relates to the peak energy of the secondary ${e}^{\pm }$ that results from either ${\pi }^{\pm }\to {e}^{\pm }$ decay or from ${\pi }^{0}\to \gamma \gamma \to {e}^{\pm }$ reactions. It will be shown later that for PKS B1424-418, up to one generation of ${e}^{\pm }$ cascade is expected, so that for both channels, the energy of ${e}^{\pm }$ is ${E}_{e}\sim 0.05{E}_{p}$. Reproducing the second hump with the observed peak frequency ${\nu }_{\mathrm{pk},2}^{\mathrm{ob}}$ from synchrotron emission of the secondary pairs requires a magnetic-field strength

$$\begin{eqnarray}&&B^{\prime} {| }_{\mathrm{pk},2}=(3\times {10}^{-2}\ {\rm{G}})\left(\displaystyle \frac{{\nu }_{\mathrm{pk},2}^{\mathrm{ob}}}{{10}^{23}\ \mathrm{Hz}}\right)\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{bulk}}}{10}\right){\left(\displaystyle \frac{{E}_{\nu }^{\mathrm{ob}}}{\mathrm{PeV}}\right)}^{-2}.\end{eqnarray} \tag{ 4 }$$

The value $B{| }_{\mathrm{pk},2}$ is clearly incompatible with $B{| }_{\mathrm{pk},1}$, which means that the proton model is not viable.

LH-psyn model: If we require a proton-synchrotron origin of the high-energy hump in the SED and a neutrino emission peaked at PeV energies, the requirement on the magnetic-field strength can be derived in a similar way as Equation (3):

$$\begin{eqnarray}&&{B}_{\mathrm{psyn}}^{\prime }=(7\times {10}^{5}\,{\rm{G}})\left(\displaystyle \frac{{\nu }_{\mathrm{pk},2}^{\mathrm{ob}}}{{10}^{23}\,\mathrm{Hz}}\right)\left(\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{bulk}}}{10}\right)\left(\displaystyle \frac{{{\rm{E}}}_{\nu }^{\mathrm{ob}}}{\mathrm{PeV}}\right).\end{eqnarray} \tag{ 5 }$$

The magnetic energy density is then of the order of ${10}^{8}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$. Compared to the photon energy density for this source, ${u}_{\mathrm{phot}}\sim 6\times {10}^{-6}{L}_{46}^{\mathrm{iso}}{R}_{18}^{-2}{{\rm{\Gamma }}}_{1}^{-4}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$, it is clearly unphysical for the jet energy budget.

LH-${\boldsymbol{\pi }}$ model: This scenario is slightly more complicated, but a few generic conditions must be met:

• 1.  
  The proton-synchrotron flux must not exceed that of synchrotron radiation of primary electrons.
• 2.  
  The Y_SSC parameter, defined as the power ratio of synchrotron-self Compton emission to synchrotron emission, must not exceed unity for this source (otherwise, the model becomes the SSC model).
• 3.  
  The observed peak frequency of the high-energy hump in the SED ${\nu }_{\mathrm{pk},2}^{\mathrm{ob}}\sim {10}^{23}$ Hz must be consistent with the characteristic energy of the secondary ${e}^{\pm }$ from pγ interactions ${E}_{e}\sim 0.05\,{E}_{p}$, with E_p determined from the neutrino energy ${E}_{\nu ,\mathrm{PeV}}^{\mathrm{ob}}$.
• 4.  
  The emission from pairs through the BH process must not overshoot the observation.

Here, the SED is approximated by four segments of power-law spectra. The equations to calculate these four constraints are 1: Equation (26), 2: Equation (27), 3: Equation (18), and 4: Equation (22) from Appendix A. All four constraints are displayed in Figure 1 for a blob radius ${R}_{\mathrm{blob}}^{\prime }={10}^{18}$ cm. When ${R}_{\mathrm{blob}}^{\prime }$ increases, the boundaries of the gray and green regions (corresponding to constraints 1 and 2, respectively) move toward the lower-left corner of the panel, whereas the regions defined by constraints 3 and 4 are independent of ${R}_{\mathrm{blob}}^{\prime }$. There is no region of overlap for all four constraints, whatever the value of ${R}_{\mathrm{blob}}^{\prime }$, which rules out the LH‐π model for PKS B1424-418.

Zoom In Zoom Out Reset image size

SSC and LH-SSC models: In the Thomson scattering regime, the frequency of the scattered photon is ${\nu }_{\mathrm{pk},2}^{\prime }\simeq {\gamma }_{e}^{\prime 2}{\nu }_{\mathrm{pk},1}^{\prime }$. Scattering proceeds in the Klein–Nishina regime for electron Lorentz factors $\gamma ^{\prime} \gtrsim {\gamma }_{\mathrm{KN}}^{\prime }$, here ${\gamma }_{\mathrm{KN}}^{\prime }h{\nu }_{\mathrm{pk},1}^{\prime }\sim {m}_{e}{c}^{2}$. Combining the two expressions, we find the following relationship among the Lorentz factors of the accelerated electrons in the comoving frame ${\gamma }_{e}^{\prime }$ and the SED parameters of PKS B1424-418:

$$\begin{eqnarray}&&{\gamma }_{e}^{\prime }=0.007\,{\gamma }_{\mathrm{KN}}^{\prime }\left(\displaystyle \frac{{\nu }_{\mathrm{pk},1}^{\mathrm{ob}}}{{10}^{14}\ \mathrm{Hz}}\right)\left(\displaystyle \frac{{\nu }_{\mathrm{pk},2}^{\mathrm{ob}}}{{10}^{23}\ \mathrm{Hz}}\right){\left(\displaystyle \frac{{\rm{\Gamma }}}{10}\right)}^{-1}\lt {\gamma }_{\mathrm{KN}}^{\prime }.\end{eqnarray} \tag{ 6 }$$

We conclude that the IC scattering of primary electrons is always Thomson scattering.

The comoving photon density can be written as

$$\begin{eqnarray}&&{u}_{\mathrm{ph}}^{\prime }={\left(\displaystyle \frac{d}{{R}_{\mathrm{blob}}}\right)}^{2}{\left(\displaystyle \frac{1+z}{{{\rm{\Gamma }}}^{2}}\right)}^{2}\displaystyle \frac{\nu {F}_{\nu }}{c}.\end{eqnarray} \tag{ 7 }$$

Here, d = 4.477 Gpc is the comoving radial distance of the source (for z = 1.522 and a flat ΛCDM universe with ${\rm{\Lambda }}=0.7$), and ${u}_{b}^{\prime }={B}^{\prime 2}/8\pi$ is the comoving magnetic-field energy density. We can express the constraints on $B^{\prime}$, ${R}_{\mathrm{blob}}^{\prime }$, and Γ in terms of observed quantities as

$$\begin{eqnarray}{R}_{\mathrm{blob}}^{\prime }{\rm{\Gamma }} & \approx & (1.4\times {10}^{19}\ \mathrm{cm}){Y}_{\mathrm{SSC}}^{-1/2}{\left(\displaystyle \frac{\nu {F}_{\nu }{| }_{\mathrm{pk},1}^{\mathrm{ob}}}{{10}^{-11}\mathrm{erg}{\mathrm{cm}}^{-2}{{\rm{s}}}^{-1}}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{{\nu }_{\mathrm{pk},1}^{\mathrm{ob}}}{{10}^{14}\mathrm{Hz}}\right)}^{-2}\left(\displaystyle \frac{{\nu }_{\mathrm{pk},2}^{\mathrm{ob}}}{{10}^{23}\ \mathrm{Hz}}\right)\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&B^{\prime} \approx (9.2\times {10}^{-3}\ {\rm{G}}){\left(\displaystyle \frac{{\nu }_{\mathrm{pk},1}^{\mathrm{ob}}}{{10}^{14}\mathrm{Hz}}\right)}^{2}{\left(\displaystyle \frac{{\nu }_{\mathrm{pk},2}^{\mathrm{ob}}}{{10}^{23}\mathrm{Hz}}\right)}^{-1}{\left(\displaystyle \frac{{\rm{\Gamma }}}{10}\right)}^{-1}\,,\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{Y}_{\mathrm{SSC}}\approx \displaystyle \frac{\nu ^{\prime} {u}_{\nu }^{\prime }({\nu }_{\mathrm{pk},1}^{\prime })}{{u}_{b}^{\prime }}\approx \displaystyle \frac{\nu {F}_{\nu }{| }_{\mathrm{pk},2}}{\nu {F}_{\nu }{| }_{\mathrm{pk},1}}\approx 10\ \end{eqnarray} \tag{ 10 }$$

and ${u}_{\nu }^{\prime }$ is the differential energy density of photons in the comoving frame. The scaling values of the observed quantities are typical for the SED of PKS B1424-418 (c.f., Figure 4). Therefore, Equation (8) implies that the blob radius is large, on the order of a light-year. The magnetic-field strength must be rather low, suggesting that the blob is located far away from the central engine, beyond the broad-line region and the dusty torus (see also Tavecchio et al. 2013). This scenario justifies the fact that we ignore the IC scattering of external photons.

The optical depth of MeV-band γ-rays to pair production, ${\tau }_{\gamma \gamma }$, can be estimated as

$$\begin{eqnarray}{\tau }_{\gamma \gamma } & \approx & \displaystyle \frac{{\rm{L}}^{\prime} }{{R}_{\mathrm{blob}}^{\prime }}\,\displaystyle \frac{{\sigma }_{T}}{4\pi {m}_{e}{c}^{3}}\approx 3\times {10}^{-5}\,\left(\displaystyle \frac{\nu {F}_{\nu }^{\mathrm{ob}}}{{10}^{-11}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{\rm{\Gamma }}}{10}\right)}^{-4}{\left(\displaystyle \frac{{R}_{\mathrm{blob}}^{\prime }}{{10}^{18}\mathrm{cm}}\right)}^{-1},\end{eqnarray} \tag{ 11 }$$

indicating that MeV γ-rays can escape the source. A numerical calculation on the optical depth shows that ${\tau }_{\gamma \gamma }$ approaches unity for multi-TeV to PeV γ-rays for typical model parameters for PKS B1424-418, which is consistent with one cascade generation of these γ-rays.

The position and flux of the high-energy peak in the SED need to be explained by the SSC model. Determining the permitted abundance of hadrons as well as the detailed fitting of the entire SED and the neutrino data require the numerical modeling of the source, which is discussed in the following section.

PKS B1424-418 exhibits significant time variabilities over a wide range of scales, as other FSRQs and BL Lacs do. Variability timescales ${t}_{\mathrm{var}}^{\mathrm{ob}}$ as short as a month are consistent with the causality argument ${t}_{\mathrm{var}}^{\mathrm{ob}}\sim {R}_{\mathrm{blob}}^{\prime }(1+z)/{\rm{\Gamma }}c$ and the constraint set by Equation (8) for a blob size ${R}_{\mathrm{blob}}=7.5\times {10}^{17}$ cm and a bulk Lorentz factor ${\rm{\Gamma }}=35$, for which ${t}_{\mathrm{var}}^{\mathrm{ob}}\simeq 3$ weeks.⁴ Faster variability has been observed on timescales shorter than a day, which may be explained by compact substructures in the jet, such as recollimation (Bromberg & Levinson 2009), "jet in a jet" scenarios, or magnetic reconnection (Giannios et al. 2009; Giannios 2013). In any case, explaining this very fast variability is beyond the scope of this paper.

We simulate time-dependent particle spectra for ${e}^{\pm }$, p, n, γ, and ${\nu }_{\alpha }$ (α denotes the neutrino flavor) by numerically solving the time-dependent differential-integral kinematic-equation system in the energy space γ, for all the particle species mentioned above:

$$\begin{eqnarray}{\partial }_{t}n(\gamma ,t) & = & -{\partial }_{\gamma }\{\dot{\gamma }(\gamma ,t)n(\gamma ,t)-{\partial }_{\gamma }[D(\gamma ,t)n(\gamma ,t)]/2\}\\ & & -\alpha (\gamma ,t)n(\gamma ,t)+Q(\gamma ,t),\end{eqnarray} \tag{ 12 }$$

where $n(\gamma )\equiv {d}^{2}N/d\gamma {dV}$ is the differential number density of the particle species. In the equation above, the source term Q may depend on energy γ, time t, and the current target-particle distributions $\{{n}_{\mathrm{tar}}(\gamma )\}$. It models the injection, emission, or generation of a new particle after an interaction, or a redistribution or re-injection of the same particle after scattering. The sink term α depends on the same set of variables and functionals, which models the escape, decay, annihilation, and disappearance of the particle from the blob due to an interaction. The differential terms $\dot{\gamma }(\gamma ,t)$ and $D(\gamma ,t)$ account for the particle cooling and diffusion effect in momentum space due to synchrotron radiation, Thomson scattering, and the BH process. In those processes, the electron or proton loses a tiny fraction of energy after scattering, and the finiteness of the numerical grid spacing is unable resolve this tiny shift in the redistribution function via the terms α and Q. Therefore, we apply the "continuous-loss" approximation, demanding an accuracy up to second-order differentiation. Due to isotropy and spherical symmetry, only the radial component $D(\gamma ,t)$ of the diffusion tensor appears in the equation.⁵

The rates and redistribution functions are described by physics, where we consider the synchrotron, IC, pair production and annihilation, photohadronic ($p\gamma$) interaction, and BH (photo-pair) processes. See Appendix B for details and Hümmer et al. (2010, hereafter H10) for the simplified $p\gamma$ interaction model, which we apply in this paper.

Bremsstrahlung and pp collisions are ignored in our case, which can be relevant in blazars when the emission region is compact (Eichmann et al. 2012). However, in our case, the large blob size leads to a low number density of cold protons, which can be easily estimated as ${n}_{p,\mathrm{cold}}\sim 2.5\times {10}^{-5}{\eta }_{p,\mathrm{cold}}{R}_{18}^{-2}{{\rm{\Gamma }}}_{1.5}^{-4}{L}_{46}^{\mathrm{ob}}$. Here, L^ob is the photon luminosity in the observer frame, and we have parameterized the energy density of cold protons as a multiple ${\eta }_{p,\mathrm{cold}}$ of photons. Adopting the zeroth-order approximation on the pp cross-section, ${\sigma }_{\mathrm{pp}}\sim 50\,\mathrm{mb}$, the optical depth of pp collision is estimated o be ${\tau }_{\mathrm{pp}}\sim {10}^{-12}{\eta }_{p,\mathrm{cold}}\,{R}_{18}^{-1}{{\rm{\Gamma }}}_{1.5}^{-4}{L}_{46}^{\mathrm{ob}}$, which is negligible. The energy-loss rate due to bremsstrahlung is estimated by Eichmann et al. (2012) to be ${\dot{\gamma }}_{\mathrm{brem}}/\gamma \sim 1.4\times {10}^{-16}(\mathrm{ln}2\,-1/3){n}_{p,\mathrm{cold}}\,\sim$ $\,1.3\,\times \,{10}^{-21}{\eta }_{p,\mathrm{cold}}{R}_{18}^{-2}{{\rm{\Gamma }}}_{1.5}^{-4}{L}_{46}^{\mathrm{ob}}$, which is significantly below that of the synchrotron loss rate ${\dot{\gamma }}_{\mathrm{syn}}/\gamma \sim {10}^{-12}\,{B}_{-3}^{2}{\gamma }_{3}$ even for the lowest-energy electrons in our case.

We use the finite-difference method to solve the equation numerically, on an evenly spaced logarithmic grid in energy and a linear one in time. The "Crank–Nicolson" differential scheme is used in time with the "Chang & Cooper" (Chang & Cooper 1970) scheme in energy to achieve stability and a more accurate goal in the correct steady-state solution. For the compatibility with the latter scheme as well as increased accuracy, we calculate up to the second-order differential term from the physics. See Appendix C for details and Vurm & Poutanen (2009, hereafter VP09) for the application to leptonic processes. Our list of input parameters and assumptions is summarized in Table 2, which we describe in greater detail in this and the next sections.

Table 2.  List of Parameters and Notations

	Group	Symbol	Definition
		${R}_{\mathrm{blob}}^{\prime }$	Comoving radius of blob, fixed to $7.5\times {10}^{17}$ cm
	Global	fesc	${e}^{\pm }$ and p escape fraction, fixed to 1/10
		${{\rm{\Gamma }}}_{\mathrm{bulk}}$	Bulk Lorentz factor of the blob fixed to 35
		$B^{\prime}$	Magnetic field strength, blob frame
		${L}_{e,\mathrm{inj}}$	Injection luminosity of primary ${e}^{-}$, AGN frame
Parameters	Leptonic	${\gamma }_{e,\min }^{\prime }$	Minimum Lorentz factor of primary ${e}^{-}$, blob frame
		${\gamma }_{e,\max }^{\prime }$	Maximum Lorentz factor of primary ${e}^{-}$, blob frame
		${\alpha }_{e,\mathrm{idx}}^{\prime }$	Power-law index of injected primary ${e}^{-}$
		${\alpha }_{p,\mathrm{idx}}^{\prime }$	Power-law index of injected protons, fixed to −2.0
	Hadronic	${\eta }_{b}$	Luminosity ratio p to ${e}^{-}$ at injection, i.e., ${\eta }_{b}\equiv {L}_{p,\mathrm{inj}}/{L}_{e,\mathrm{inj}}$
		${E}_{p,\max }^{\mathrm{ob}}$	Maximum energy of injected protons, observer frame
	${P}_{\mathrm{0,1,0}}({E}_{p,\max },{\eta }_{b})$		Probability to observe 0, 1, and 0 neutrino events in the 0.5–1.6 PeV, 1.6–2.4 PeV, and >2.4 PeV bands in IceCube, respectively, as a function of ${E}_{p,\max }$ and ${\eta }_{b}$
	${P}_{\nu ,\max }$		Maximum value of ${P}_{\mathrm{0,1,0}}({E}_{p,\max },{\eta }_{b})$ in the parameter space, named "neutrino best fit"
	${L}_{\nu }$		Total neutrino luminosity, including all flavors
	${L}_{\gamma }$		γ-ray luminosity integrated over the frequency band ${10}^{18.1}\sim {10}^{24.4}$ Hz
Notations	+		SED best-fit mark
	×		Neutrino best-fit mark, global
	⨂		Neutrino best fit under the constraint of the SED being reproduced within $3\sigma$ confidence
	*		Joint best-fit mark, for SED and neutrino
	(1)		Neutrino flux to expect one $\gt 0.5$ PeV event during the 2LAC phase, assuming the same effective area of IceCube as in the IC-2yr and Burst phases (see Figure 5)

Download table as:  ASCIITypeset image

The dynamical SEDs of PKS B1424-418 reported in K16 are categorized into four phases: (1) a flare in 2010, lasting about one month, (2) 2LAC phase, from 2008.8 to 2010.9, (3) IC-2yr period, from 2010.5 to 2012.5, the first two years of IceCube observations; (4) Burst phase, from 2012.6 to 2013.3, when the source experienced a long-lasting high-flux phase in γ-rays. Figure 2 provides a visual timeline, including the time-averaged GeV-band γ-ray flux. A 2 PeV neutrino event in IceCube (IC35, also dubbed "Big Bird") was observed on 2012 December 4, during phase 4, with a position consistent with that of PKS B1424-418. The SEDs shown in K16 are based on time-averaged spectra from each phase. In this paper, the flare phase is ignored since its duration was too short to result in a neutrino fluence comparable with that of other phases.⁶ We therefore focus on three phases (2LAC, IC-2yr, and Burst), as indicated in Figure 2, which can be interpreted as the time-dependent evolution of the AGN. We will first simulate the phases independently, and then interpret the evolution (changes) of the parameters.

Zoom In Zoom Out Reset image size

For each phase, the SED and neutrino spectra are modeled by a steady-state solution to Equation (12). Under the SSC and LH-SSC scenarios, the first peak of the SED is described as synchrotron emission from primary ${e}^{-}$, and the γ-rays are described as SSC emission from the same ${e}^{-}$ population (see Table 1). After fixing a few global parameters such as R_blob and ${{\rm{\Gamma }}}_{\mathrm{bulk}}$, the following two-step simulations are performed for each of the three phases of PKS B1424-418 (2LAC, IC-2yr, and Burst). (1) Use leptonic simulations $({\eta }_{b}\equiv 0)$ to find the best-fit parameters of the primary ${e}^{-}$ for the low-energy hump and γ-ray band (${10}^{22}\mbox{--}{10}^{25}$ Hz). (2) Inject protons until their spectrum reaches a steady state to find the total SED and neutrino spectrum.

In step 1, with leptonic simulations, the following parameter space is scanned: ${L}_{e,\mathrm{inj}}({10}^{42.5}\sim {10}^{45.5}$ erg s⁻¹) $\otimes \,{\gamma }_{e,\min }^{\prime }({10}^{2.6}\sim {10}^{3.9})$ $\otimes \,{\gamma }_{e,\max }^{\prime }({10}^{4.2}\sim {10}^{6.0})$ $\,\otimes \,{\alpha }_{e,\mathrm{inj}}(-2.0\sim -1.0)$ $\,\otimes \,B^{\prime} ({10}^{-2.2}\,\sim {10}^{-4.0}\,{\rm{G}})$ $\otimes \,{{\rm{\Gamma }}}_{\mathrm{bulk}}(1\sim 200)$. The best-fit parameters are obtained by ${\chi }^{2}$-minimization.

In step 2, including protons, the SED and the neutrino spectrum are calculated for a logarithmic ${\eta }_{b}({10}^{2}\sim {10}^{8})\,\otimes {E}_{p,\max }^{\prime }({10}^{2}\sim {10}^{9}\,\mathrm{GeV})$ parameter grid for each phase of PKS B1424-418. Here, ${\eta }_{b}$ is the proton-to-electron luminosity ratio (baryonic loading) at injection, ${\eta }_{b}={L}_{p,\mathrm{inj}}/{L}_{e,\mathrm{inj}}$, and ${E}_{p,\max }^{\prime }$ is the maximum proton energy. Each point of the grid corresponds to a unique combination of ${\eta }_{b}$ and ${E}_{p,\max }$, and therefore an independent simulation. The entire grid has a resolution of 400 × 160 for each phase, requiring a total of 192,000 hadronic simulations.

The global parameters ${R}_{\mathrm{blob}}^{\prime }=7.5\times {10}^{17}$ cm and ${{\rm{\Gamma }}}_{\mathrm{bulk}}=35$ are chosen according to Equation (8), including a variability time ${t}_{\mathrm{var}}\lt 1$ month. The energy-independent escape rate for ${e}^{-}$ and p is assumed to be ${t}_{e,p}^{-1}=0.1\,{t}_{\mathrm{fs}}^{-1}$, which is an order of magnitude slower than the free-streaming escape rate (that of the neutrinos). We also use a fixed power-law injection index ${\alpha }_{p,\mathrm{inj}}=-2.0$ for proton injection, where the results are not very sensitive to this value.

TANAMI VLBI data of PKS B1424-418 indicate an angular size of a few milli-arcsecond for the emission region in the radio band, which translates into a physical size of about the order of a few 10¹⁹ cm, which is larger than that of the blob. The extended VLBI component can be interpreted as synchrotron emission from electrons that escaped from the blob into an extended region of size ${R}_{\mathrm{ext}}\sim 3.0\times {10}^{19}$ cm that is filled with a weaker magnetic field ${B}_{\mathrm{ext}}^{\prime }\simeq 0.1\,{B}_{\mathrm{blob}}^{\prime }$. We assume that the radio data are fitted by synchrotron emission from electrons that escaped the blob into an extended region using these parameters. All other contributions to the SED come from the blob.

In this subsection, we focus on the Burst phase and the relationship to the potentially observed neutrino event during that phase.

For each parameter-space simulation, we perform independent optimizations. First, we adapt our model to the SED and find the best fit and confidence regions by calculating the reduced ${\chi }^{2}$ values.⁷ This SED best fit is marked with the symbol "+." Then, with the predicted neutrino spectrum from each simulation, we calculate the probability, ${P}_{010}({E}_{p,\max },{\eta }_{b})$ of observing 0, 1, and 0 neutrino events from PKS B1424-418 with IceCube, as indeed observed in the (0.5–1.6), (1.6–2.4), and $\gt 2.4$ PeV energy bands, respectively. In each band, the expected number of neutrino events follows a Poisson distribution. The best adaptation to the neutrino data without regard for the SED is marked by the symbol "×," and its fit probability is denoted by ${P}_{\nu ,\max }$. The joint best-fit point, by maximizing the joint probability of the SED and neutrino fits, is marked with the symbol "✳."

We find that both SED and neutrino fits independently prefer large baryonic loadings for the Burst phase, which may point toward a baryonically loaded burst. At the neutrino and joint best-fit points in Figure 3 (×, ✳), the associated proton maximum energies are around 10 PeV, which is consistent with the observed PeV neutrino event. However, the maximum proton flux at the neutrino best fit, ×, is in tension with the X-ray data, which we will demonstrate below.

Zoom In Zoom Out Reset image size

We show the SED for the SED best fit (+), the joint best fit (✳), and the neutrino best fit (×) in the left panel of Figure 3. We clearly observe that the SED is in tension with the data in the X-ray energy range. On the other hand, the SED is described reasonably well for the other two fit points. Note that the radio data are fitted from the extended emission region.

As a next step, we address the question whether the observed neutrino event can come from the Burst phase of PKS B1424-418, as reported in K16. One test of this hypothesis is energetics—in K16, the neutrino luminosity was directly related to the section of the SED we defined "${L}_{\gamma }$" in the caption of Figure 4, One can easily see that this energy range is dominated by leptonic processes (1) in our model, contrary to what has been assumed in K16 (for which hadronic processes 2–4 would need to dominate that energy range).

Zoom In Zoom Out Reset image size

We list the fractional contributions of the different physical interactions to the SED for the joint best-fit case of the Burst phase in the caption of Figure 4, where one can clearly read off that the SSC contribution dominates. We also show the neutrino-to-γ-ray ratio ${L}_{\nu }/{L}_{\gamma }$, which is of order 5% for the models fitting the SED. These numbers are to be interpreted as an additional "theory" correction factor in addition to those included in K16 (independent of the spectral correction in K16, which would reduce this number). They reflect the fact that hadronic processes only dominate in a small portion of the energy range marked "${L}_{\gamma }$" in Figure 4, considering that the second hump cannot be dominated by hadronic processes. Note that hadronic components contribute significantly to the SED outside the pre-defined energy range of ${L}_{\gamma }$ (${10}^{18.1}\sim {10}^{24.4}\,\mathrm{Hz}$, 5 keV to 10 GeV). In K16, the predicted number of neutrino events in the 1.0–2.0 PeV bin is 1.6, assuming that the entire second hump is generated from hadronic processes. From our model, this needs to be corrected by this factor of 0.05, arriving at ∼0.08 events. Consistently, the neutrino spectrum from our numerical simulation (the joint best-fit case) predicts 0.094 events in IceCube within the same energy bin.

Let us now address if we can draw a self-consistent picture of the AGN blazar over time. For that purpose, we independently fit the parameters for the three phases 2LAC, IC-2yr, and Burst in Figure 2; these fits are shown in Figures 3, 5, and 6. For the IC-2yr phase, we predict neutrino events in the same format as the one in the Burst phase, even though no neutrinos are detected during this phase. The purpose is to demonstrate that the correlation between the PeV neutrino event with the Burst phase is weak. Indeed, the IC-2yr phase predicts up to a probability ${P}_{\mathrm{0,1,0}}=5.7 \%$ of producing the same observations in IceCube, compared to ${P}_{\mathrm{0,1,0}}=3.2 \%$ in the Burst phase (see Table 3). The relative probability is therefore ${P}_{\mathrm{0,1,0}}(\mathrm{IC} \mbox{-} 2\mathrm{yr}):{P}_{\mathrm{0,1,0}}(\mathrm{Burst})={\rm{1.6:1}}$. Here ${P}_{\mathrm{0,1,0}}$ is defined as the joint probability of observing 0, 1, and 0 neutrino events in the 0.5–1.6, 1.6–2.4, and $\gt 2.4$ PeV energy bins, respectively, where in each bin, the expected number of neutrino events follows a Poisson distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Parameters and the Expected Number of Neutrino Events in IceCube for the Best-fit Models During Each Phase

Phase	2LAC			IC-2yr			Burst
Time	2008.9–2010.9			2010.5–2012.5			2012.6–2013.3
${R}_{\mathrm{blob}}/\mathrm{cm}$	$7.5\times {10}^{17}$
${{\rm{\Gamma }}}_{\mathrm{bulk}}$	35
$B^{\prime} /\mathrm{mG}$	2.5			2.0			2.5
${L}_{e,\mathrm{inj}}/{L}_{\mathrm{edd}}$	$6.7\times {10}^{-5}$			$1.2\times {10}^{-4}$			$3.2\times {10}^{-4}$
${\gamma }_{e,\min }^{\prime }$	$1.8\times {10}^{3}$			$1.8\times {10}^{3}$			$2.2\times {10}^{3}$
${\gamma }_{e,\max }^{\prime }$	$1.3\times {10}^{5}$			$8.9\times {10}^{4}$			$1.0\times {10}^{5}$
${\alpha }_{e,\mathrm{idx}}^{\prime }$	−2.2			−2.2			−1.8
${\alpha }_{p,\mathrm{idx}}^{\prime }$	−2.0
Fit symbol	+	⨂	(1)	+	*	×	+	*	×
Fit	SED	Table 2	Table 2	SED	joint	ν	SED	joint	ν
${\eta }_{b}/{10}^{6}$	11.5	1.8	2.2	2.3	0.76	0.72	0.33	0.10	0.20
${E}_{p,\max }^{\mathrm{ob}}/\mathrm{PeV}$	0.37	7.5	8.3	0.68	6.8	11.2	0.68	7.5	12.3
${N}_{\nu },0.5\mbox{--}1.6\,\mathrm{PeV}$	0	0.35	0.74	0.0090	0.62	0.79	0.0045	0.21	0.76
${N}_{\nu },1.6\mbox{--}2.4\,\mathrm{PeV}$	0	0.082	0.18	0	0.13	0.22	0	0.042	0.22
${N}_{\nu },\,\,\,\,\gt 2.4\,\mathrm{PeV}$	0	0.037	0.10	0	0.044	0.21	0	0.018	0.24
${P}_{\mathrm{0,1,0}}$, %	⋯	⋯	⋯	0	5.7	6.5	0	3.2	6.4
Photon SED ${\chi }^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.$	0.83	1.36	86.9	1.82	1.94	7.10	1.38	1.40	171

Note. ${P}_{\mathrm{0,1,0}}$ next to the bottom row is the joint probability of observing 0, 1, and 0 neutrino events in the three energy bins above. See Table 2 for the definition of symbols and notations.

Download table as:  ASCIITypeset image

For the 2LAC phase, IceCube was operating with mainly the IC40+59 configuration, compared to the full IC86 for the other phases. Here we show the contours of the predicted neutrino events of $\gt 0.5$ PeV during this phase in the parameter space, from the IC40+59 configuration. The best-fit point for the SED in Figure 5 is marked with "+," the maximum neutrino flux within the 3σ region for the SED fit is marked by "⨂," and a representative point on the ${N}_{\nu }=1$ contour is picked and marked with "①." Since the SED of the 2LAC phase has a lower photon flux than that of the IC-2yr or Burst phase, a lower target photon density in the source is implied. In order to have a sufficiently large $p\gamma$ interaction rate to produce the right amount of X-rays and one neutrino event, the required proton-to-lepton luminosity ratio ${\eta }_{b}$ has to be very large: ${\eta }_{b}\gtrsim {10}^{6}$. From that perspective, it is not surprising that no neutrinos were observed during that phase.

The X-ray band, seen as the gap between the two humps in the SED, is a mixture of leptonic and hadronic components, as shown in Figure 4. Hadronic secondary emission depends on the detailed processing of the electromagnetic cascades in the source, mainly via the following channels: $p\gamma \,\to {\pi }^{0}\to \gamma \gamma \to {e}^{\pm }$, $p\gamma \to {\pi }^{\pm }\to {\mu }^{\pm }\to {e}^{\pm }$, and $p\gamma \to p+{e}^{\pm }$. In the ${\pi }^{0}$ channel, the emission predominantly comes from the first generation of pairs. The source itself becomes optically thick for γ-rays above roughly 100 TeV. It was further estimated by Stecker et al. (2012) that the spectrum above ∼100 GeV will suffer from additional suppression due to absorption by EBL, but this effect does not play a significant role in our fitting.

We list the detailed best-fit parameters and expected number of neutrino events during each phase in Table 3. We note that the parameters for the magnetic fields, and the minimum and maximum energies for the ${e}^{-}$ and p injections, are very similar. The SEDs evolving from the earliest phase, 2LAC, to IC-2yr, and finally to the Burst phase can be simply reproduced by an increase of the ${e}^{-}$ versus $p$ injection rate with time, and a slight hardening effect on the ${e}^{-}$ injection spectrum.

The proton-to-lepton luminosity ratio required to fit the neutrino result is large, of the order of 10⁵–10⁶. This value is coincidentally similar to that found for other blazars, such as by Diltz et al. (2015) and Petropoulou et al. (2015). In our case, on one hand, the pγ interaction rate is lower than those in the above cases and more protons are needed to achieve the same level of $\gamma$-rays; on the other hand, compared to their LH-$\pi$ models, we only need a subdominant hadronic contribution to $\gamma$-rays, which lowers the requirement for proton luminosity. The physical proton-injection luminosities needed for these values are around a few times the Eddington luminosity. This requirement can be alleviated if we assume a lower escape rate for protons, e.g., a magnetic configuration that can trap the protons longer in the blob. In our model, due to the low $p\gamma$ efficiency, most protons simply escape without having a $p\gamma$ interaction and the interactions cause no visible effects on the proton spectrum. Therefore, the required ${L}_{p,\mathrm{inj}}$ scales linearly with the escape rate.

If we consider the relativistic electrons only, i.e., ${\gamma }_{e}\gt 10$, the differential cross-section (expressed in units of ${\sigma }_{T}$) can be greatly simplified under the head-on collision approximation (Dermer & Menon 2009, hereafter DM09):

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\sigma }_{\mathrm{IC}}}{d{\gamma }_{f,o}d\mu }({\gamma }_{e,o},{\gamma }_{e,i},{\gamma }_{f,i},\mu )=\displaystyle \frac{3}{8}\displaystyle \frac{1}{{\gamma }_{e,i}{x}^{3}}\left[{{yx}}^{2}+1+2x+\displaystyle \frac{1}{y}({x}^{2}-2x-2)+\displaystyle \frac{1}{{y}^{2}}\right]\,H\left(y-\displaystyle \frac{1}{1+2x}\right)H\left(1-\displaystyle \frac{x}{2{\gamma }_{e,i}^{2}}-y\right),\end{eqnarray} \tag{ 34 }$$

where $x={\gamma }_{e,i}{\gamma }_{f,i}(1-\mu )$ is the photon energy in the electron rest frame and $y=1-{\gamma }_{f,o}/{\gamma }_{e,i}$. The subindex i(o) stands for the incoming (outgoing) particle and e(f) represents electron (photon). ${\gamma }_{e}$ and ${\gamma }_{f}$ are dimensionless energy units defined by ${\gamma }_{e}={E}_{e}/{m}_{e}{c}^{2}$ and ${\gamma }_{f}={E}_{f}/{m}_{e}{c}^{2}$,

$$\begin{eqnarray}&&{R}_{\mathrm{IC}}({\gamma }_{f,o}\leftarrow {\gamma }_{e,i},{\gamma }_{f,i})=\displaystyle \frac{c}{2}\displaystyle \int (1-\mu )\displaystyle \frac{{d}^{2}{\sigma }_{\mathrm{IC}}}{d{\gamma }_{f,o}d\mu }({\gamma }_{e,o},{\gamma }_{e,i},{\gamma }_{f,i},\mu )d\mu \\ =\,\left\{\begin{array}{ll}\displaystyle \frac{6(1-u)\tfrac{u}{v}\mathrm{ln}\left[\tfrac{u}{2v(1-u)}\right]-3\left(\tfrac{u}{v}-2+2u\right)\left(\tfrac{u}{v}+1+\tfrac{{u}^{2}}{2}-u\right)}{4v{\gamma }_{e,i}{(u-1)}^{3}}, & \,\mathrm{for}\,\displaystyle \frac{u}{2{\gamma }_{e,i}^{2}}\lt v\lt \displaystyle \frac{2u}{1+2u}\,\mathrm{and}\,u\lt 2{\gamma }_{e,i}^{2}-1/2\\ 0 & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 35 }$$

where $u=2{\gamma }_{e,i}$ and $v={\gamma }_{f,o}/{\gamma }_{e,i}$.

In the Thomson scattering limit, it reduces to

$$\begin{eqnarray}&&{R}_{\mathrm{IC}}({\gamma }_{f,o}\leftarrow {\gamma }_{e,i},{\gamma }_{f,i})=\displaystyle \frac{\mathrm{ln}(w/2)-(w/2-1/2-{w}^{-1})}{2{\gamma }_{f,o}/3},\end{eqnarray} \tag{ 36 }$$

where $w=u/v=2{\gamma }_{e,i}^{2}{\gamma }_{f,i}/{\gamma }_{f,o}$.

The total reaction rate, averaged over the angle μ, is

$$\begin{eqnarray}&&{R}_{\mathrm{IC}}({\gamma }_{f,i},{\gamma }_{e,i})=\displaystyle \int {R}_{\mathrm{IC}}({\gamma }_{f,o}\leftarrow {\gamma }_{e,i},{\gamma }_{f,i})d{\gamma }_{f,o}\\ &&\qquad \,\qquad =\,\displaystyle \frac{3c}{8{u}^{3}(1+2u)}[-2u(4+9u+{u}^{2})+(4+u){(1+2u)}^{2}\mathrm{ln}(1+2u)+4u(1+2u)\mathrm{Li}(2,-2u)],\end{eqnarray} \tag{ 37 }$$

which has the asymptotic forms

$$\begin{eqnarray}{R}_{\mathrm{IC}}({\gamma }_{f,i},{\gamma }_{e,i})=\,\left\{\begin{array}{ll}c\left(1-\displaystyle \frac{4u}{3}\right), & \mathrm{for}\,u\ll 1(\mathrm{Thomson})\\ c\displaystyle \frac{3}{4u}\left[\mathrm{ln}\left(\displaystyle \frac{u}{2}\right)-\displaystyle \frac{1}{2}\right], & \mathrm{for}\ u\gg 1(\mathrm{Klein}\mbox{--}\mathrm{Nishina})\end{array}\right..\end{eqnarray} \tag{ 38 }$$

$\mathrm{Li}(2,x)$ is the dilogarithm defined as

$$\begin{eqnarray}&&\mathrm{Li}(2,z)=-{\int }_{0}^{z}\displaystyle \frac{\mathrm{ln}(1-t)}{t}{dt}.\end{eqnarray} \tag{ 39 }$$

In the Thomson regime, after scattering, the electron loses a tiny fraction of energy, and therefore the function ${R}_{\mathrm{IC}}({\gamma }_{e,o}\leftarrow {\gamma }_{f,i}{\gamma }_{e,i})$ is highly peaked around ${\gamma }_{e,i}$, which the numerical grid is unable to resolve. Here, we use the differential terms to account for this effect in the numerical computation, up to the second order (VP09):

$$\begin{eqnarray}&&-{\partial }_{\gamma }\{{\dot{\gamma }}_{\mathrm{IC}}n({\gamma }_{e},t)-{\partial }_{\gamma }[{D}_{\mathrm{IC}}({\gamma }_{e},t)n({\gamma }_{e},t)]\}=-{\alpha }_{\mathrm{IC}}({\gamma }_{e},t)n({\gamma }_{e},t)+{Q}_{\mathrm{IC}}({\gamma }_{e},t).\end{eqnarray} \tag{ 40 }$$

Using the moment expansion (Equations (C11), (C12), (C18), and (C19) of VP09), they are expressed as

$$\begin{eqnarray}&&{\dot{\gamma }}_{\mathrm{IC},e}({\gamma }_{e})=\int {\gamma }_{f}({{\rm{\Psi }}}_{1}-{{\rm{\Psi }}}_{0}){n}_{f}({\gamma }_{f})d{\gamma }_{f}\,,\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{D}_{\mathrm{IC},e}({\gamma }_{e})=\int {\gamma }_{f}^{2}({{\rm{\Psi }}}_{2}-2{{\rm{\Psi }}}_{1}+{{\rm{\Psi }}}_{0}){n}_{f}({\gamma }_{f})d{\gamma }_{f},\end{eqnarray} \tag{ 42 }$$

where these moments are given by Nagirner & Poutanen (1994),

$$\begin{eqnarray}{{\rm{\Psi }}}_{0}({\gamma }_{f},{\gamma }_{e}) & \approx & 1-\displaystyle \frac{2}{3}(4{\gamma }_{e}^{2}-1){\gamma }_{f}{\gamma }_{e}^{-1}+\displaystyle \frac{26}{5}(2{\gamma }_{e}^{2}-1){\gamma }_{f}^{2}\\ {{\rm{\Psi }}}_{1}({\gamma }_{f},{\gamma }_{e}) & \approx & \displaystyle \frac{1}{3}(4{\gamma }_{e}^{2}-1)-\displaystyle \frac{1}{5}(42{\gamma }_{e}^{4}-29{\gamma }_{e}^{2}+2){\gamma }_{f}{\gamma }_{e}^{-1}+\displaystyle \frac{1}{25}(1176{\gamma }_{e}^{4}-1147{\gamma }_{e}^{2}+206){\gamma }_{f}^{2}\\ {{\rm{\Psi }}}_{2}({\gamma }_{f},{\gamma }_{e}) & \approx & \displaystyle \frac{1}{15}(42{\gamma }_{e}^{4}-34{\gamma }_{e}^{2}+7)-\displaystyle \frac{4}{75}(528{\gamma }_{e}^{6}-618{\gamma }_{e}^{4}+172{\gamma }_{e}^{2}-7){\gamma }_{f}{\gamma }_{e}^{-1}\\ & & +\displaystyle \frac{1}{525}(109120{\gamma }_{e}^{6}-158856{\gamma }_{e}^{4}+63677{\gamma }_{e}^{2}\,-\,6066){\gamma }_{f}^{2}.\end{eqnarray} \tag{ 43 }$$

For the pair production process, ${\gamma }_{f,1}+{\gamma }_{f,2}\to {\gamma }_{e,1}+{\gamma }_{e,2}$, two photons with dimensionless energies ${\gamma }_{f,i}={E}_{f,i}/{m}_{e}{c}^{2},(i=1,2)$ are annihilated and an electron−positron pair is produced with Lorentz factors ${\gamma }_{e,1},{\gamma }_{e,2}$, respectively. Obviously, for energy conservation, we have ${\gamma }_{f,1}\,+{\gamma }_{f,2}={\gamma }_{e,1}+{\gamma }_{e,2}$. We use the treatment and expressions from VP09:

$$\begin{eqnarray}{Q}_{\mathrm{pp}}({\gamma }_{e,1}) & = & c{\displaystyle \int }_{{\gamma }_{f,1}^{* }}^{\infty }n({\gamma }_{f,1})d{\gamma }_{f,1}{\displaystyle \int }_{{\gamma }_{f,2}^{* }}^{\infty }{R}_{\mathrm{pp}}({\gamma }_{e,1}\leftarrow {\gamma }_{f,1},{\gamma }_{f,2})n({\gamma }_{f,2})d{\gamma }_{f,2}\end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray}{R}_{\mathrm{pp}}({\gamma }_{e,1}\leftarrow {\gamma }_{f,1},{\gamma }_{f,2}) & = & -{\gamma }_{f,1}^{-2}{\gamma }_{f,2}^{-2}[S({\gamma }_{e,1},{\gamma }_{f,2},{\gamma }_{f,1},{w}_{U})\,-S({\gamma }_{e,1},{\gamma }_{f,2},{\gamma }_{f,1},{w}_{L})]/4,\end{eqnarray} \tag{ 45 }$$

in which

$$\begin{eqnarray}&&S({\gamma }_{e,1},{\gamma }_{f,2},{\gamma }_{f,1},w)=-{[{({\gamma }_{f,1}+{\gamma }_{f,2})}^{2}-4{w}^{2}]}^{1/2}+T({\gamma }_{e,1},{\gamma }_{f,2},{\gamma }_{f,1},w)+T({\gamma }_{e,1},{\gamma }_{f,1},{\gamma }_{f,2},w),\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}&&T(x,y,z,w)={w}^{3}{({yz})}^{-3/2}({yz}-1){h}^{-1}[{A}_{0}(h)-\,{(1+h)}^{1/2}]-{(1+h)}^{1/2}{w}^{-1}{({yz})}^{-1/2}\\ &&\qquad \qquad \,\qquad +\,(w/2){({yz})}^{-3/2}[{(1+h)}^{-1/2}({y}^{2}+{yz}+{xz}-{xy}-2{w}^{2})\,-\,4{{yzA}}_{0}(h)],\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}{A}_{0}(h)\left\{\begin{array}{ll}=\,{h}^{-1/2}\mathrm{ln}[{h}^{1/2}+{(1+h)}^{1/2}] & \mathrm{for}\ h\gt 0\\ =\,{(-h)}^{-1/2}\arcsin \sqrt{-h} & \mathrm{for}\ h\lt 0\\ \approx \,1-h/6+3{h}^{2}/40 & \mathrm{for}\ h\approx 0\end{array}\right.,\end{eqnarray} \tag{ 48 }$$

and finally, with $h=[{(x-y)}^{2}-1]{w}^{2}/{yz}$ and the integration boundaries ${w}_{L}={w}_{-}$, ${w}_{U}=\min [\sqrt{{\gamma }_{f,1}{\gamma }_{f,2}},{w}_{+}]$, in which

$$\begin{eqnarray}&&{w}_{\pm }=[{\gamma }_{e,1}{\gamma }_{e,2}+1\pm \sqrt{({\gamma }_{e,1}^{2}+1)({\gamma }_{e,2}^{2}+1)}]/2.\end{eqnarray} \tag{ 49 }$$

The emergence rate for photons, due to the pair annihilation process, is

$$\begin{eqnarray}{Q}_{\mathrm{pa}}({\gamma }_{f,1}) & = & c{\displaystyle \int }_{{\gamma }_{e,2}^{* }}^{\infty }n({\gamma }_{e,2})d{\gamma }_{e,2}{\displaystyle \int }_{{\gamma }_{e,1}^{* }}^{\infty }n({\gamma }_{e,1}){R}_{\mathrm{pa}}({\gamma }_{f,1}\leftarrow {\gamma }_{e,1},{\gamma }_{e,2}).\end{eqnarray} \tag{ 50 }$$

The expression for R_pa can be obtained via the symmetry

$$\begin{eqnarray}{R}_{\mathrm{pa}}({\gamma }_{f,1}\leftarrow {\gamma }_{e,1},{\gamma }_{e,2}) & = & {R}_{\mathrm{pa}}({\gamma }_{f,2}\leftarrow {\gamma }_{e,1},{\gamma }_{e,2})\,=\,{R}_{\mathrm{pp}}({\gamma }_{e,1}\leftarrow {\gamma }_{f,1},{\gamma }_{f,2})\,=\, & {R}_{\mathrm{pp}}({\gamma }_{e,2}\leftarrow {\gamma }_{f,1},{\gamma }_{f,2})\end{eqnarray} \tag{ 51 }$$

and energy conservation

$$\begin{eqnarray}&&{\gamma }_{f,1}+{\gamma }_{f,2}={\gamma }_{e,1}+{\gamma }_{e,2}.\end{eqnarray} \tag{ 52 }$$

The lower limits of the integration in Equation (44) are

$$\begin{eqnarray}{\gamma }_{f,1}^{* } & = & \displaystyle \frac{1}{2}{\gamma }_{e,1}(1-{\beta }_{e,1})\\ {\gamma }_{f,2}^{* } & = & \left\{\begin{array}{ll}{\gamma }_{f,2}/\{[2{\gamma }_{f,2}-{\gamma }_{e,1}(1+{\beta }_{e,1})]{\gamma }_{e,1}(1+{\beta }_{e,1})\} & \,\mathrm{for}\,x\gt {x}_{+}\\ {\gamma }_{f,2}/\{[2{\gamma }_{f,2}-{\gamma }_{e,1}(1-{\beta }_{e,1})]{\gamma }_{e,1}(1-{\beta }_{e,1})\} & \,\mathrm{for}\,x\lt {x}_{-}\\ {\gamma }_{e,1}-{\gamma }_{f,2}+1 & \,\mathrm{for}\,{x}_{-}\leqslant x\leqslant {x}_{+}\end{array}\right.\end{eqnarray} \tag{ 53 }$$

with ${x}_{\pm }=[1+{\gamma }_{e,1}(1\pm {\beta }_{e,1})]/2$, while the limits in Equation (50) are

$$\begin{eqnarray}{\gamma }_{e,2}^{* } & = & \left\{\begin{array}{ll}{\gamma }_{A} & \,\mathrm{for}\,{\gamma }_{f,2}\lt 1/2\\ 1 & \,\mathrm{for}\,{\gamma }_{f,2}\geqslant 1/2\end{array}\right.{\,\gamma }_{e,1}^{* } & = & \left\{\begin{array}{ll}{\gamma }_{-} & \mathrm{for}\ {\gamma }_{f,2}\leqslant 1/2\\ {\gamma }_{-} & \mathrm{for}\ 1/2\lt {\gamma }_{f,2}\lt 1\,\mathrm{and}\,{\gamma }_{+}\lt {\gamma }_{B}\\ {\gamma }_{+} & \mathrm{for}\ {\gamma }_{f,2}\geqslant 1\,\mathrm{and}\,{\gamma }_{+}\lt {\gamma }_{B}\\ 1 & \,\mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 54 }$$

where

$$\begin{eqnarray}{\gamma }_{\pm } & = & ({F}_{\pm }+{F}_{\pm }^{-1})/2,\,{F}_{\pm }=2{\gamma }_{f,2}-{\gamma }_{e,2}(1\pm {\beta }_{e,2}),\,{\gamma }_{A} & = & {\gamma }_{f,2}+1/(4{\gamma }_{f,2}),{\gamma }_{B}={\gamma }_{f,2}-({\gamma }_{f,2}-1)/(2{\gamma }_{f,2}-1).\end{eqnarray} \tag{ 55 }$$

The disappearance rate for photons ${\gamma }_{f,2}$ due to pair production with a target photon ${\gamma }_{f,1}$ is

$$\begin{eqnarray}&&{\alpha }_{\mathrm{pp}}({\gamma }_{f,2})=\int {R}_{\mathrm{pp}}({\gamma }_{f,2},{\gamma }_{f,1})n({\gamma }_{f,1})d{\gamma }_{f,1},\end{eqnarray} \tag{ 56 }$$

and R_pp (in units of ${\sigma }_{T}$) is given by (see also DM09)

$$\begin{eqnarray}&&R({\gamma }_{f,2},{\gamma }_{f,1})=\displaystyle \frac{3}{8}c{\gamma }_{f,2}^{-2}{\gamma }_{f,1}^{-2}\bar{\varphi }(u,v)\end{eqnarray} \tag{ 57 }$$

and

$$\begin{eqnarray}&&\bar{\varphi }(u,v)=[2v+{(1+v)}^{-1}]\mathrm{ln}u-{\mathrm{ln}}^{2}u-2(2v+1){v}^{1/2}\,{(v+1)}^{-1/2}+4\mathrm{ln}u\mathrm{ln}(1+u)+{\pi }^{2}/3+4{\mathrm{Li}}_{2}(-u),\end{eqnarray} \tag{ 58 }$$

where

$$\begin{eqnarray}&&v={\gamma }_{f,2}{\gamma }_{f,1}-1,u=\displaystyle \frac{\sqrt{v+1}+\sqrt{v}}{\sqrt{v+1}-\sqrt{v}}.\end{eqnarray} \tag{ 59 }$$

The asymptotic form of Equation (58) is

$$\begin{eqnarray}\bar{\varphi }(u,v)\approx \left\{\begin{array}{ll}2v\mathrm{ln}(4v)-4v+{\mathrm{ln}}^{2}(4v) & \mathrm{for}\ v\gg 1\\ \displaystyle \frac{4}{3}{v}^{3/2} & \mathrm{for}\ v\ll 1\end{array}.\right.\end{eqnarray} \tag{ 60 }$$

The pair annihilation rate is

$$\begin{eqnarray}&&{\alpha }_{\mathrm{pa}}({\gamma }_{e,2})=\int {R}_{\mathrm{pa}}({\gamma }_{e,1},{\gamma }_{e,2})n({\gamma }_{e,1})d{\gamma }_{e,1},\end{eqnarray} \tag{ 61 }$$

where the kernel, in units of ${\sigma }_{T}$, is

$$\begin{eqnarray}&&{R}_{\mathrm{pa}}({\gamma }_{e,1},{\gamma }_{e,2})=\displaystyle \frac{3}{8}c{\gamma }_{e,1}^{-2}{\gamma }_{e,2}^{-2}[{S}_{\mathrm{pa}}({\gamma }_{e}^{+})-{S}_{\mathrm{pa}}({\gamma }_{e}^{-})],\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{S}_{\mathrm{pa}}(x)={x}^{2}\mathrm{ln}(4{x}^{2})-2{x}^{2}+\displaystyle \frac{3}{4}{\mathrm{ln}}^{2}(4{x}^{2}),\end{eqnarray} \tag{ 63 }$$

and

$$\begin{eqnarray}&&{\gamma }_{e}^{\pm }=\displaystyle \frac{1}{2}[{\gamma }_{e,1}{\gamma }_{e,2}+1\pm \sqrt{({\gamma }_{e,1}^{2}-1)({\gamma }_{e,2}^{2}-1)}].\end{eqnarray} \tag{ 64 }$$

The emission rate of synchrotron photons by a relativistic electron ${\gamma }_{e}\gtrsim 10$ is given by

$$\begin{eqnarray}&&{Q}_{\mathrm{syn}}({\gamma }_{f})=\int {R}_{\mathrm{syn}}({\gamma }_{f}\leftarrow {\gamma }_{e})n({\gamma }_{e})d{\gamma }_{e}.\end{eqnarray} \tag{ 65 }$$

The integration kernel is well known as

$$\begin{eqnarray}&&{R}_{\mathrm{syn}}({\gamma }_{f}\leftarrow {\gamma }_{e})=c{\sigma }_{T}(3\sqrt{3}/\pi ){\gamma }_{f}^{-1}{u}_{b}{F}_{b}^{-1}{z}^{2}\,\{{K}_{4/3}(z){K}_{1/3}(z)-(3/5)z[{K}_{4/3}^{2}(z)-{K}_{1/3}^{2}(z)]\},\end{eqnarray} \tag{ 66 }$$

where

$$\begin{eqnarray}{u}_{b}={B}^{2}/(8\pi {m}_{e}{c}^{2}),{F}_{b}=B/{B}_{\mathrm{crit}},{B}_{\mathrm{crit}} & = & {({m}_{e}{c}^{2})}^{2}/({ce}{\hslash })=4.41\times {10}^{13}\,{\rm{G}},z= & {\gamma }_{f}{\gamma }_{e}^{-2}{F}_{b}^{-1}/3,\end{eqnarray} \tag{ 67 }$$

and K_n(z) is the modified Bessel function of the second kind.

For the extinction rate of the photons due to the synchrotron-self-absorption effect, we follow the treatment of VP09:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{ssa}}({\gamma }_{f})=\displaystyle \frac{{\lambda }_{C}^{3}}{8\pi }\int {\gamma }_{f}^{-1}{R}_{\mathrm{syn}}({\gamma }_{f}\leftarrow {\gamma }_{e}){\gamma }_{e}^{2}{\partial }_{\gamma }[{\gamma }_{e}^{-2}n({\gamma }_{e})]d{\gamma }_{e},\end{eqnarray} \tag{ 68 }$$

where ${\lambda }_{C}=h/{m}_{e}c$ is the Compton wavelength. The cooling effect on electrons, positrons, and protons, due to synchrotron and synchrotron-self-absorption effects, is modeled as a continuous energy-loss process and thus described by the differential terms

$$\begin{eqnarray}{\dot{\gamma }}_{\mathrm{syn}}(\gamma ) & = & {\dot{\gamma }}_{s}(\gamma )+2{\gamma }^{-1}{H}_{s}(\gamma )+{\partial }_{\gamma }{H}_{s}(\gamma ),{D}_{\mathrm{syn}}(\gamma ) & = & 2{H}_{s}(\gamma ),\end{eqnarray} \tag{ 69 }$$

and

$$\begin{eqnarray}&&{\dot{\gamma }}_{s}=-\displaystyle \frac{4}{3}{\sigma }_{T}{u}_{B}{\gamma }_{e}^{2}\,,\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{H}_{s}({\gamma }_{e})=\displaystyle \frac{{\lambda }_{C}^{3}}{8\pi }\int {R}_{\mathrm{syn}}({\gamma }_{f}\leftarrow {\gamma }_{e})n({\gamma }_{f})d{\gamma }_{f}.\end{eqnarray} \tag{ 71 }$$

For the $p\gamma$ interaction $p+f\to X$, with p being a proton here, f the target photon, and X representing a proton, neutron, or pion, we follow the simplified treatment sim A in H10, and incorporate them in this numerical framework. The generation rate of particle X can be expressed as

$$\begin{eqnarray}&&{Q}_{X}({E}_{X})=\int {n}_{f}({E}_{f}){R}_{p\gamma }({E}_{X}\leftarrow {E}_{f}){{dE}}_{f},\end{eqnarray} \tag{ 72 }$$

where $Q\equiv d{\dot{n}}_{X}/{{dE}}_{X}$ and $n\equiv {d}^{2}N/{dVdE}$ so that the phenomenological values in the references can be directly applied.

The response function ${R}_{p\gamma }({E}_{X}\leftarrow {E}_{f})$ depends on a further layer of integration:

$$\begin{eqnarray}&&{R}_{p\gamma }({E}_{X})=\int {R}_{p\gamma }({E}_{X}\leftarrow {E}_{p},{E}_{f}){n}_{p}({E}_{p}){{dE}}_{p},\end{eqnarray} \tag{ 73 }$$

where the integration kernel is

$$\begin{eqnarray}&&{R}_{p\gamma }({E}_{X}\leftarrow {E}_{p},{E}_{f})={{cE}}_{p}^{-1}\displaystyle \sum _{i}\delta (x-{\chi }_{i}){M}_{i}{f}_{i}(y),\end{eqnarray} \tag{ 74 }$$

where i stands for the interaction channel, M_i is the multiplicity of the particle X, ${\chi }_{i}$ is the inelasticity of the collision, f(y) is defined as an integration over the cross-section as a function of photon energy ${\epsilon }_{r}$ expressed in the particle rest frame,

$$\begin{eqnarray}&&f(y)=\displaystyle \frac{1}{2{y}^{2}}{\int }_{{\epsilon }_{{th}}}^{2y}d{\epsilon }_{r}{\epsilon }_{r}\sigma ({\epsilon }_{r}),\end{eqnarray} \tag{ 75 }$$

and finally, with $x={E}_{X}/{E}_{p}$ and $y={\gamma }_{p}{E}_{f}$.

Under the $\delta$-function approximation, ${R}_{p\gamma }({E}_{X}\leftarrow {E}_{f})$ simplifies to

$$\begin{eqnarray}&&{R}_{p\gamma }({E}_{X}\leftarrow {E}_{f})=c\displaystyle \sum _{i}{E}_{p,0}^{-1}{\chi }_{i}^{-2}{E}_{X}{M}_{i}{f}_{i}({y}_{0}){n}_{p}({E}_{p,0}),\end{eqnarray} \tag{ 76 }$$

where ${E}_{p,0}$ and y₀ corresponds to the solution of $x-{\chi }_{i}=0$. The expression for Q_X simplifies to

$$\begin{eqnarray}&&{Q}_{X}({E}_{X})={{cn}}_{p}({\chi }_{i}^{-1}{E}_{X}){E}_{X}^{-1}{m}_{p}\displaystyle \sum _{i}{M}_{i}\int {n}_{f}({E}_{X}^{-1}{\chi }_{i}{m}_{p}y){f}_{i}(y){dy}.\end{eqnarray} \tag{ 77 }$$

The disappearance rate, due to participation in the $p\gamma$ process for protons, is

$$\begin{eqnarray}{\alpha }_{p\gamma }({E}_{p})\,=\,\displaystyle \int {{dE}}_{X}{R}_{p\gamma }({E}_{X}\leftarrow {E}_{p})\,=\,\displaystyle \int {{dE}}_{X}\displaystyle \int {{dE}}_{f}{R}_{p\gamma }({E}_{X}\leftarrow {E}_{p},{E}_{f}){n}_{f}({E}_{f})\,=\, & c\displaystyle \sum _{i}{M}_{i}\displaystyle \int {{dE}}_{f}{n}_{f}({E}_{f}){f}_{i}({\gamma }_{p}{E}_{f}),\end{eqnarray} \tag{ 78 }$$

and for photons,

$$\begin{eqnarray}{\alpha }_{p\gamma }({E}_{f}) & = & \displaystyle \int {{dE}}_{X}{R}_{p\gamma }({E}_{X}\leftarrow {E}_{f})\,= & c\displaystyle \sum _{i}{\chi }_{i}^{-1}{M}_{i}\displaystyle \int {{dE}}_{X}{f}_{i}({\chi }_{i}^{-1}{E}_{X}{E}_{f}{m}_{p}^{-1}){n}_{p}({\chi }_{i}^{-1}{E}_{X}).\end{eqnarray} \tag{ 79 }$$

The expressions of f_i(y) and coefficients are given by Equations (30), (33)–(35), and (40) and Tables 3, 5, 6 of HU10 when i falls into the categories of resonance, direct, and multi-pion production, respectively.

For the BH or photo-pair process $p+{\gamma }_{f}\to {\gamma }_{e,1}+{\gamma }_{e,2}$, we incorporate the calculation using the framework similar to the multi-pion production in the previous section. The generation rate for pairs is given by Equation (77), where the multiplicity M = 1 for ${e}^{-}$ and ${e}^{+}$, respectively. The cross-section and inelasticity are approximated by a series of step functions, where each step is defined as an interaction channel i. The cross-section ${\sigma }_{\mathrm{BH}}({\gamma }_{f}^{\prime })$ as a function of photon energy (in units of ${m}_{e}{c}^{2}$) ${\gamma }_{f}^{\prime }$ in the proton rest frame, obtained under the approximation of zero recoil of the proton, can be expressed as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{BH},\mathrm{tot}}({\gamma }_{f}^{\prime })={\int }_{1}^{{\gamma }_{f}^{\prime }}d{\gamma }_{e,1}{\int }_{-1}^{+1}d{\mu }_{e,1}\displaystyle \frac{{d}^{2}{\sigma }_{\mathrm{BH},\mathrm{diff}}}{d{\gamma }_{e,1}d{\mu }_{e,1}},\end{eqnarray} \tag{ 80 }$$

where the differential cross-section in the proton rest frame is expressed as (Blumenthal 1970):

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\sigma }_{\mathrm{BH},\mathrm{diff}}}{d{\gamma }_{e,1}d{\mu }_{e,1}}=\left(\displaystyle \frac{3{\alpha }_{f}{\sigma }_{T}{p}_{1}{p}_{2}}{16\pi {\gamma }_{f}^{3}}\right)\left[-4(1-{\mu }^{2})\displaystyle \frac{2{\gamma }_{1}^{2}+1}{{p}_{1}^{2}{{\rm{\Delta }}}_{1}^{4}}\right.+\displaystyle \frac{5{\gamma }_{1}^{2}-2{\gamma }_{1}{\gamma }_{2}+3}{{p}_{1}^{2}{{\rm{\Delta }}}_{1}^{2}}+\displaystyle \frac{{p}_{1}^{2}-{\gamma }_{f}^{2}}{{T}^{2}{{\rm{\Delta }}}_{1}^{2}}+\displaystyle \frac{2{\gamma }_{2}}{{p}_{1}^{2}{{\rm{\Delta }}}_{1}}+\displaystyle \frac{Y}{{p}_{1}{p}_{2}}\phantom{\rule{}{3.65ex}}(2{\gamma }_{1}(1-{\mu }^{2})\\ &&\,\times \left.\displaystyle \frac{3{\gamma }_{f}+{p}_{1}^{2}{\gamma }_{2}}{{{\rm{\Delta }}}_{1}^{4}}+\displaystyle \frac{2{\gamma }_{1}^{2}({\gamma }_{1}^{2}+{\gamma }_{2}^{2})-7{\gamma }_{1}^{2}-3{\gamma }_{1}{\gamma }_{2}-{\gamma }_{2}^{2}+1}{{{\rm{\Delta }}}_{1}^{2}}+\displaystyle \frac{{\gamma }_{f}({\gamma }_{1}^{2}-{\gamma }_{1}{\gamma }_{2}-1)}{{{\rm{\Delta }}}_{1}}\right)\\ &&\,\left.-\displaystyle \frac{{\delta }_{+}^{T}}{{p}_{2}T}\left(\displaystyle \frac{2}{{{\rm{\Delta }}}_{1}^{2}}-\displaystyle \frac{3{\gamma }_{f}}{{{\rm{\Delta }}}_{1}}-\displaystyle \frac{{\gamma }_{f}({p}_{1}^{2}-{\gamma }_{f}^{2})}{{T}^{2}{{\rm{\Delta }}}_{1}}\right)-\displaystyle \frac{2{y}_{+}}{{{\rm{\Delta }}}_{1}}\right],\end{eqnarray} \tag{ 81 }$$

where

$$\begin{eqnarray}&&T=| {\boldsymbol{k}}-{{\boldsymbol{p}}}_{1}| ,\,Y=\displaystyle \frac{2}{{p}_{1}^{2}}\mathrm{ln}\left[\displaystyle \frac{{\gamma }_{1}{\gamma }_{2}+{p}_{1}{p}_{2}+1}{{\gamma }_{f}}\right],{y}_{+}={p}_{2}^{-1}\mathrm{ln}\left[\displaystyle \frac{{\gamma }_{2}+{p}_{2}}{{\gamma }_{2}-{p}_{2}}\right],\,{\delta }_{\,+}^{T}=\mathrm{ln}\left[\displaystyle \frac{T+{p}_{2}}{T-{p}_{2}}\right],{{\rm{\Delta }}}_{1}={\gamma }_{1}(1-{\beta }_{1}\mu ),\end{eqnarray} \tag{ 82 }$$

and in the above equations we have dropped the primes and the electron subindex e (in which, $1={e}^{-},2={e}^{+}$).

The inelasticity for one electron can be obtained via

$$\begin{eqnarray}&&{\chi }_{\mathrm{BH}}({\gamma }_{f})=\displaystyle \frac{{m}_{e}}{{m}_{p}}{\int }_{1}^{{\gamma }_{f}-1}d{\gamma }_{1}{\int }_{-1}^{+1}d\mu {{\rm{\Delta }}}_{1}\displaystyle \frac{{d}^{2}{\sigma }_{\mathrm{BH},\mathrm{diff}}}{d{\gamma }_{1}d\mu }.\end{eqnarray} \tag{ 83 }$$

Equations (80) and (83) are approximated as a sum of step functions:

$$\begin{eqnarray}{\sigma }_{\mathrm{BH},\mathrm{tot}}({\gamma }_{f}) & = & \displaystyle \sum _{i}{\sigma }_{\mathrm{BH},i,0}H({\gamma }_{f}-{\gamma }_{f,i}^{\min })H({\gamma }_{f,i}^{\max }-{\gamma }_{f})\\ {\chi }_{\mathrm{BH}}({\gamma }_{f}) & = & \displaystyle \sum _{i}{\chi }_{\mathrm{BH},i,0}H({\gamma }_{f}-{\gamma }_{f,i}^{\min })H({\gamma }_{f,i}^{\max }-{\gamma }_{f}),\end{eqnarray} \tag{ 84 }$$

where

$$\begin{eqnarray}{\sigma }_{\mathrm{BH},i,0}({\gamma }_{f}) & = & \displaystyle \frac{1}{{\gamma }_{f,i}^{\max }-{\gamma }_{f,i}^{\min }}{\displaystyle \int }_{{\gamma }_{f,i}^{\min }}^{{\gamma }_{f,i}^{\max }}d{\gamma }_{f}{\sigma }_{\mathrm{BH},\mathrm{tot}}({\gamma }_{f})\\ {\chi }_{\mathrm{BH},i,0}({\gamma }_{f}) & = & \displaystyle \frac{1}{{\gamma }_{f,i}^{\max }-{\gamma }_{f,i}^{\min }}{\displaystyle \int }_{{\gamma }_{f,i}^{\min }}^{{\gamma }_{f,i}^{\max }}d{\gamma }_{f}{\chi }_{\mathrm{BH}}({\gamma }_{f}),\end{eqnarray} \tag{ 85 }$$

and the functions {${f}_{\mathrm{BH},i}(y)$} are obtained from Equation (40) of HU10 by replacing the coefficients with the ones from the above equation.

The extinction function of the photons is

$$\begin{eqnarray}{\alpha }_{\mathrm{BH}}({\gamma }_{f}) & = & 2c\displaystyle \sum _{i}{\chi }_{\mathrm{BH},i,0}^{-1}\displaystyle \int {{dE}}_{1}{f}_{\mathrm{BH},i}({\chi }_{\mathrm{BH},i,0}^{-1}{E}_{1}{E}_{f}{m}_{p}^{-1}){n}_{p}({\chi }_{\mathrm{BH},i,0}^{-1}{E}_{1}).\end{eqnarray} \tag{ 86 }$$

For protons, the generation rate in principle can be obtained from Equation (77) by substituting ${M}_{i}=1$, ${f}_{i}={f}_{\mathrm{BH},i}$, and ${\chi }_{i}=1-2{\chi }_{\mathrm{BH},i,0}$. However, since ${\chi }_{\mathrm{BH},i,0}\ll 1$, ${\chi }_{i}\approx 1$, the term ${n}_{p}({\chi }_{i}^{-1}{E}_{p})$ almost overlaps with the parent proton bin ${n}_{p}({E}_{p})$, and the numerical grid cannot resolve this difference. Instead, similar to the case of IC scattering in the Thomson regime, this effect is treated as a continuous energy-loss process and thus described by the differential terms. By requiring equality between the integral terms and the differential terms, such as Equation (40), we have

$$\begin{eqnarray}{\dot{E}}_{\mathrm{BH}}({E}_{p}) & = & \displaystyle \sum _{i}-2{\chi }_{\mathrm{BH},i,0}{E}_{p}{\alpha }_{p,i}({E}_{p})\\ {D}_{\mathrm{BH}}({E}_{p}) & = & \displaystyle \sum _{i}4{\chi }_{\mathrm{BH}}^{2}{E}_{p}^{2}{\alpha }_{p,i}({E}_{p}),\end{eqnarray} \tag{ 87 }$$

where we have defined ${\alpha }_{p,i}({E}_{p})=2c\int {{dE}}_{f}{n}_{f}({E}_{f}){f}_{\mathrm{BH},i}({\gamma }_{p}{E}_{f})$.

The dominant decay channels of ${\pi }^{\pm }$ and ${\mu }^{\pm }$ are

$$\begin{eqnarray}{\pi }^{+}({\pi }^{-}) & \to & {\mu }^{+}({\mu }^{-})+{\nu }_{\mu }({\bar{\nu }}_{\mu })\\ {\mu }^{+}({\mu }^{-}) & \to & {e}^{+}({e}^{-})+{\nu }_{e}({\bar{\nu }}_{e})+{\bar{\nu }}_{\mu }({\nu }_{\mu }).\end{eqnarray} \tag{ 88 }$$

The probability density function for the decay product j from the parent particle i, as a function of scaling variable $x={E}_{j}/{E}_{i}$, is (Lipari et al. 2007)

$$\begin{eqnarray}{f}_{{\mu }_{R}^{+}}(x) & = & {f}_{{\mu }_{L}^{-}}(x)=\displaystyle \frac{{r}_{\mu \pi }^{2}(1-x)}{{(1-{r}_{\mu \pi }^{2})}^{2}x}H(x-{r}_{\mu \pi }^{2})\\ {f}_{{\mu }_{L}^{+}}(x) & = & {f}_{{\mu }_{R}^{-}}(x)=\displaystyle \frac{x-{r}_{\mu \pi }^{2}}{{(1-{r}_{\mu \pi }^{2})}^{2}x}H(x-{r}_{\mu \pi }^{2}).\end{eqnarray} \tag{ 89 }$$

Due to energy conservation, the corresponding neutrino directly from pion decay has a distribution function of

$$\begin{eqnarray}&&{f}_{\nu }(1-x)={f}_{\mu }(x).\end{eqnarray} \tag{ 90 }$$

For ${\pi }^{0}\to \gamma +\gamma$, the distribution function for a photon is

$$\begin{eqnarray}&&{f}_{\gamma }(x)=\displaystyle \frac{1}{2{\beta }_{\pi }{\gamma }_{\pi }}H\left(x-\displaystyle \frac{1-{\beta }_{\pi }}{2}\right)H\left(\displaystyle \frac{1+{\beta }_{\pi }}{2}-x\right).\end{eqnarray} \tag{ 91 }$$

The neutrino distribution function from muon decay is

$$\begin{eqnarray}{f}_{{\bar{\nu }}_{\mu }}(x,h) & = & {f}_{{\nu }_{\mu }}(x,-h)=\left(\displaystyle \frac{5}{3}-3{x}^{2}+\displaystyle \frac{4}{3}{x}^{3}\right)+h\times \left(-\displaystyle \frac{1}{3}+3{x}^{2}-\displaystyle \frac{8}{3}{x}^{3}\right)\\ {f}_{{\nu }_{e}}(x,h) & = & {f}_{{\bar{\nu }}_{e}}(x,-h)=(2-6{x}^{2}+4{x}^{3})+h\times (2-12x+18{x}^{2}-8{x}^{3}),\end{eqnarray} \tag{ 92 }$$

and for ${e}^{\pm }$, since we no longer need to distinguish between their chiralities in this paper, the distribution can be simply expressed as (e.g., Particle Data Group)

$$\begin{eqnarray}&&{f}_{e}(x)=\displaystyle \frac{4}{3}(1-{x}^{3})\end{eqnarray} \tag{ 93 }$$

under the relativistic approximation, ${\gamma }_{e}\gtrsim 10$.