A closer look at the electromagnetic signatures of Bethe-Heitler pair production process in blazars

The “twin birth” of a positron and an electron by a photon in the presence of a nucleus, known as Bethe-Heitler pair production, is a key process in astroparticle physics. The Bethe-Heitler process offers a way of channeling energy stored in a population of relativistic protons (or nuclei) to relativistic pairs with extended distributions. Contrary to accelerated leptons, whose maximum energy is limited by radiative losses, the maximal energy of pairs is determined by the kinematics of the process and can be as high as the parent proton energy. We take a closer look at the features of the injected pair distribution, and provide a novel empirical function that describes the spectrum of pairs produced by interactions of single-energy protons with single-energy photons. The function is the kernel of the Bethe-Heitler pair production spectrum that replaces a double numerical integration involving the complex differential cross section of the process, and can be easily implemented in numerical codes. We further examine under which conditions Bethe-Heitler pairs produced in blazar jets can emit γ-ray photons via synchrotron radiation, thus providing an alternative to the inverse Compton scattering process for high-energy emission in jetted active galactic nuclei. For this purpose, we create 36 numerical models using the code ATHEνA optimized so that the Bethe-Heitler synchrotron emission dominates their γ-ray emission. After taking into consideration the broadband spectral characteristics of the source, the jet energetics, and the properties of radiation fields present in the blazar environment, we conclude that γ-rays in high-synchrotron-peaked blazars are unlikely to be produced by Bethe-Heitler pairs, because the emitting region is found to be opaque in photon-photon pair production at photon energies ≳ 10 GeV. On the contrary, γ-ray spectra of low-synchrotron-peaked blazars may arise from Bethe-Heitler pairs in regions of the jet with typical transverse size ∼ 1015 – 1016 cm and co-moving magnetic field 50 – 500 G. For such cases, an external thermal target photon field with temperatures T ∼ 4 · 102– 6 · 103 K is needed. The latter values could point to the dusty torus of the AGN. Interestingly, a Bethe-Heitler-dominated high-energy component is mostly found in models of intermediate-synchrotron peaked blazars, for a wide range of magnetic fields and source radii.


Introduction
The "twin birth" of a positron and an electron by a photon in the presence of a nucleus, known as Bethe-Heitler pair production, is a key process in astroparticle physics.While the cross section of the interaction has been provided by Hans Bethe and Walter Heitler in 1934 [1], it was not until the discovery of the cosmic-wave background radiation in the mid-sixties [2] and the detection of ultrahigh-energy cosmic rays [3], that Bethe-Heitler pair production received attention by the astrophysics community, and was discussed as an energy loss mechanism for cosmic rays during their propagation to Earth [e.g.[4][5][6].
With the discovery of X-ray non-thermal compact sources, such as active galactic nuclei (AGN), there was a focus shift to pair cascades [e.g.7,8].Because acceleration of electrons to relativistic energies can be radiation-limited, other processes that can efficiently channel energy to relativistic pairs have been discussed, including Bethe-Heitler pair production [e.g.[9][10][11].Because this is a subdominant energy loss channel for relativistic protons that satisfy the energy threshold for pion production, Bethe-Heitler pair injection was often neglected in AGN radiation models.Nevertheless, Bethe-Heitler pairs are injected with a different energy distribution than the pairs from photopion production interactions [e.g.12].Therefore, Bethe-Heitler pairs may still leave a radiative signature on the AGN spectrum as illustrated in more recent numerical works [13][14][15][16][17].
To the best of our knowledge no analytical description of the pair production spectrum arising from the interaction of a single proton with a single photon exists in the literature.Numerical codes for multi-messenger emission that take into account the injection of pairs from Bethe-Heitler interactions [e.g.[18][19][20][21][22] employ numerical integration of at least two integrals 1 , one involving the energyand angle-dependent cross section of the process [6].In works aiming at the study of non-linear pair cascades (where the target photon field evolves with time and depends on the parent proton distribution) or in cases where high-accuracy is needed, the calculation of the Bethe-Heitler spectrum can be a computational bottleneck (for ways to speed up the computation, see Refs.[21] and [22]).
In this work we construct an empirical analytical function for the Bethe-Heitler pair injection spectrum that can be easily implemented in numerical codes (section 2).Our function is benchmarked with numerical results from the code ATHEνA that are based on Monte-Carlo simulations of proton-photon interactions performed by [23].The function is not a unique mathematical representation of the pair production spectrum.Instead, it was conceptualized by studying the shape of the pair production spectrum, and identifying scaling relations with the proton Lorentz factor and the target photon energy as seen in the proton rest frame (i.e. the interaction energy).Using the empirical function we take a closer look at key features of the pair injection spectrum produced in interactions between protons and photons with power-law distributions, which are typically expected in astrophysical sources.By considering these key features of the injection spectrum, we examine under which conditions Bethe-Heitler pairs can produce γ-ray spectra as those observed from AGN jets (section 3), and conclude with a discussion of our findings (section 4).

The Bethe-Heitler pair production spectrum
The distribution of electrons or positrons produced by a single proton with Lorentz factor γ p interacting with a population of photons with distribution f ph (ϵ) is given by [12]: where W(ω, E − , ξ) is the cross section of the interaction [6], E − is the electron energy, ω is the photon energy in the proton rest frame, and ξ ≡ γ p E − −γ e γ p p − where γ p and γ e are the proton and pair Lorentz factors, respectively.In the expression above2 all energies are normalized to m e c 2 .
In codes that model the time evolution of particle distributions the pair distribution at injection has to be computed at every time step [18,19,24,25].Numerical integration of the last integral involving the cross section can slow down the computations, especially if these aim at high accuracy.The last two integrals can be instead pre-calculated for various combinations of γ e , γ p , and ω values [19,25].By using interpolation in a three-dimensional space, the desired integral values can then be obtained by reducing the computational time by a factor of ∼ 20 − 30 [21].Nevertheless, the computation would be even faster, if an empirical analytical function describing the pair distribution was available.A similar approach was adopted by [12], who derived analytical parametrizations for the energy distributions of secondary particles produced in proton-photon interactions that lead to the production of pions.
Our goal is to provide an analytical function for the Bethe-Heitler pair production spectrum (section 2.1) that could be easily implemented in numerical codes or used in semi-analytical calculations.Moreover, by using this empirical function we will study key features of the pair production spectrum, such as characteristic energies and power-law slopes (section 2.2) for various energy distributions of interacting particles.

An empirical function for the pair production spectrum
We present an analytical function that describes the total injection spectrum of particles (i.e., electrons and positrons, later referred to as pairs), q BH (γ e ), produced by the Bethe-Heitler interaction process of a single proton of Lorentz factor γ p and a single target photon of energy ϵ (in units of m e c 2 ).Our function is benchmarked with numerical results from ATHEνA [18,26] that are based on Monte-Carlo simulations of proton-photon interactions performed by [23].Details about the function design and its performance against ATHEνA results for different cases are presented in Appendix A. We also compare the results of our novel function, q BH , against those expected by the numerical implementation of the Bethe-Heitler production spectrum presented in [12] (see Appendix A.3).
Our goal is to describe the injection spectrum of secondaries produced by the interaction of a single proton with a single photon.In ATHEνA, where all particle energies are discretized, monoenergetic distributions are approximated by particles occupying a single bin of energy.More specifically, ATHEνA employs logarithmic energy grids with resolution of 0.1 for protons and pairs, and 0.2 for photons.For example, a mono-energetic proton distribution with Lorentz factor 10 5 in ATHEνA is described by a power-law distribution of slope s p in the logarithmic bin [5.0, 5.1].Furthermore, a single mono-energetic proton in ATHEνA is described by the proton distribution in the selected energy bin divided by the total number of protons occupying this bin.
Since ATHEνA is used to benchmark the empirical Bethe-Heitler injection spectrum, the latter function is bound to inherit the limitations of its numerical implementation.For instance, the Monte Carlo simulations performed by [23], and used in tabulated form in ATHEνA, account for interaction energies up to 10 4 (i.e., γ p ϵ ≤ 10 4 ).Furthermore, the high-energy tail of the pair production spectrum (γ e ≳ γ p ) has saw-like edges due to low-number statistics of the Monte Carlo simulations (e.g., see figure 1c).The aforementioned edges make it difficult to model the spectrum and introduce uncertainties in the description of the high-energy cutoff.Finally, we benchmark the empirical functions for interactions involving protons with γ p > 10.The reason is that the self-similarity of Bethe-Heitler injection spectra for near-threshold interactions breaks down for protons with γ p < 10.All the above shape the validity regime of our empirical function q BH (γ e ).
We first present pair production spectra computed with ATHEνA (figure 1) in order to highlight the main features that we aimed at capturing with the analytical function.We observe that the differential production rate peaks at a Lorentz factor γ e,pk ≈ ϵ −1 , with ϵ being the target photon energy (see panel a in figure 1).Meanwhile, the position of the aforementioned injection rate peak is not affected by the proton Lorentz factor (see figure 18 in Appendix A).However, this behavior changes when we look at the energy distribution of the produced pairs, γ 2 e q BH (γ e ).When the interaction takes place near the threshold (i.e.γ p ϵ ≳ 2), most of energy is transferred to pairs with γ e ≈ γ p (see panel b).Both the peak Lorentz factor and the peak of the energy injection rate scale linearly with γ p for interactions close to the threshold (see also panel d).On the contrary, the pair energy distribution from interactions taking place away from the threshold tends to spread over a wider range of pair Lorentz factors without having a well-defined peak (panel c).Instead, these interactions seem to transfer roughly equal amounts of energy into pairs that are produced with a wider distribution of Lorentz factors (see redder curve in panel c).For interactions happening away from the energy threshold, the energy injection spectra show large fluctuations at high Lorentz factors due to low event statistics of the Monte Carlo simulations that are provided, in a tabular form, as input to the code.As a result, the high-energy tail of the analytical spectrum will be more uncertain (see also Appendix A).
The empirical function we constructed reads, where γ e,pk = (1.23ϵ)−1 is the Lorentz factor where the number of injected pairs peaks (i.e.where γ e q BH (γ e ) has a maximum), a 1 , a 2 , a 3 and γ e,cr are summarized in table 1 and a log 10 (γ p ϵ) (2.4) m e γ p .Branches (2), ( 3) and ( 4) of table 1 control the sharpness of the high-energy cutoff of the injection function, q BH .
The normalization of the injection rate A(γ p , ϵ) is given by (in units of s −1 ) + H log 10 (γ p ϵ) h 1 e h 2 log 10 (γ p ϵ) + J log 10 (ϵ) − 3 log 10 (R 15 ) where the parameters are listed in table 2 and R ′ b,15 is the radius of the spherical blob normalized to 10 15 cm.Additionally, p(γ p ϵ) is the power-law index appearing in the first term of the exponential Comparison of the analytical function against ATHEνA for all the bench-marking cases we used can also be found in Appendix A. We present next the performance of the analytical function against numerical results for a more realistic case that involves interactions between protons and photons with extended power-law distributions.The pair production spectrum in this case is a convolution of q BH with the two particle distributions, with N p (γ p ) = dN p / dγ p and N ph (ϵ) = dN ph / dϵ being the differential number distributions of protons and photons respectively.We consider a power-law proton distribution N p (γ p ) ∝ γ −2 p extending from γ p,min = 10 4 to γ p,max = 10 8 interacting with a power-law photon distribution, N ph (ϵ) ∝ ϵ −s ph .We adopt five different photon indices, s ph = −1, 0, 1, 2, 3, that cover a wide range of potentially relevant astrophysical sources.For example, non-thermal photon spectra from high-energy astrophysical sources are usually soft with photon indices s ph ≥ 2, while the Rayleigh-Jeans part of a black-body spectrum has s ph = −1.The photon spectrum extends from ϵ min = 10 −8 to ϵ max = 10 −4 .With the adopted energy ranges the highest energy protons will "see" as targets the whole power-law photon distribution, while the lowest energy protons will be able to interact only with the highest energy photons of the distribution.In all cases, the target photon compactness, defined as is fixed and equal to 0.1.Similarly, the proton compactness is ℓ p = 10 −6 and R ′ b = 10 15 cm.Our results are presented in figure 2. The kernel of the Bethe-Heitler pair production spectrum, q BH , when convolved with the particle power-law distributions yields pair spectra that are in good agreement with the numerical results of ATHEνA (see left panel).We observe that the peak position and the peak value of the differential pair injection rate are well reproduced.The shape of the curve for a wide range of pair Lorentz factors is overall well described for all values of the photon index we explored.Only at the high-energy tail of the distribution we start seeing a discrepancy between the numerical results and semi-analytical spectra.The maximum difference we find is about factor of 2.5 (see ratio plot in the left panel of 2), but the production rate at these Lorentz factors is many orders of magnitude lower than its peak value.These differences at high energies become more pronounced when comparing the energy distributions of produced pairs -see right panel of figure 2. For soft target photon distributions (s ph = 3), the semi-analytical result overestimates by a factor of almost 3, at most, the energy injection rate at the high-energy tail of the distribution.On the contrary, it slightly underestimates the energy injection rate by up to a similar factor when the target photon spectra are very hard (s ph ≤ 0).Nonetheless, the analytical function can reproduce the shape of the energy distribution function for a wide range of Lorentz factors, including the peak position and its peak value.

Key features of the pair production spectrum and synchrotron emission
The empirical function of Eq. (2.2) describes well the general features of the numerically computed pair production spectrum.In this paragraph, we therefore use the empirical function to study key features of the pair injection spectrum produced by interacting power-law (PL) distributions (section 2.2.1), as well as the resulting synchrotron spectra that are a common observable in non-thermal astrophysical sources (section 2.2.2).

Injection spectra of Bethe-Heitler pairs
In figure 2 (left panel) one may notice that the differential injection rate peaks at γ e,pk ≈ 1/ϵ max ≈ γ p,min regardless of the photon index, with the exception of the s ph = 3 case where the pair injection spectrum appears almost flat for γ e ∈ [ϵ −1 max , ϵ −1 min ].On the contrary, the peak position of the energy injection spectra depends on the photon index as shown in the right panel of figure 2.
Here, we explore the dependence of the energy injection spectrum on both s p and s ph by considering cases with s p = −1, 0, 1, 2, 3 and s ph = −1, 0, 1, 2. In all examples, we use ℓ ph = 10 −1 , ℓ p = 10 −6 , and R ′ b = 10 15 cm.In figure 3 we present a decomposition of the energy injection spectrum into separate components produced by protons of specific Lorentz factors.For proton distributions with s p ≤ 1 (where the number of protons, γ p dN p dγ p ∝ γ −s p +1 p , peaks at the highest energies for s p < 1 or protons are equally distributed in number across all energies (s p = 1)), we notice that the pair energy distribution is determined by the most energetic protons regardless of s ph .However, for a fixed proton slope, changes in the photon index affect the shape of the produced distribution (see e.g.second row from the top).For hard photon spectra (e.g.s ph = −1, 0) the spectrum appears to be almost flat.The reason for this behavior is that most photons are concentrated in the high-energy tail of the distribution, thus providing far-from threshold targets for the higher energy protons, which are also more numerous.The aforementioned conclusion also agrees with the shape of the pair energy distribution (for s ph = −1, 0 and s p = −1, 0, 1) which coincides with the away-from-threshold spectral shape (redder curves); see also spectra for γ p ϵ ≫ 2 in figure 1.Similarly, for softer photon spectra (i.e.larger values of s ph ) most photons are found in the low-energy part of the distribution.As a result, near-threshold interactions of the high-energy protons become in this case more important than those taking place away from threshold.The latter leads to a harder spectrum with a well defined maximum (see e.g.panel for s p = −1 and s ph = 2).
Figure 3 also shows that when low-energy protons dominate in number (i.e.s p > 1) they can determine the maximum of the pair energy distribution for certain choices of the photon index: when target photons are equally distributed across the energy spectrum (s ph = 1) or concentrated near the high-energy tail of the photon distribution (s ph = −1, 0).In these cases threshold interactions of the low-energy protons are the ones being dominant (bluer lines).Furthermore, we notice that pair production by more energetic protons is also important in those scenarios, since they are the ones creating the power-law segment of the distribution beyond its peak.Nonetheless, for a fixed value of the proton index (e.g.s p = 2 or s p = 3), we observe that as the photon index increases the contribution of the high-energy protons (redder lines) becomes more significant, thus dominating the distribution and defining its peak (cases with s p = 2 and s ph = 2) or balancing the low-energy proton contributions (cases with s p = 3 and s ph = 2).It is noteworthy that the contributions of the highenergy and low-energy protons to the total pair energy spectrum are comparable when s p − s ph ∼ 1.

Synchrotron spectra from Bethe-Heitler pairs
Here we use the analytical approximation of the Bethe-Heitler injection rate to compute the photon energy distribution emitted by the pairs.After all, the photon spectrum is the main observable of an astrophysical source.Given that synchrotron radiation is relevant to the majority of non-thermal emitting sources, we will compare the synchrotron spectrum emitted by Bethe-Heitler pairs as obtained from ATHEνA and by our empirical function.To do so, we calculate the steady-state distribution of pairs, as described in [27]: dγ ′ e Q(γ ′ e ) e γ br γ ′ e (2.8) with τ esc ∼ R ′ b /c being the escape timescale of pairs from the source and Q is computed from Eq. (2.6).Furthermore, γ br represents the cooling break Lorentz factor, which in the case of synchrotron radiative losses, is defined as with u ′ B = B ′2 /(8π) being the energy density of the magnetic field.Finally, we compute the synchrotron spectrum arising from the steady state distribution of Eq. (2.8) using the full expression for the synchrotron emissivity [28].
The results presented in figure 4 are obtained for the same parameters as those used for figure 2 and magnetic field strength B ′ = 40 G; additional cases for different magnetic field strengths are presented in Appendix B. The distribution of cooled pairs peaks (in number) at γ e ≈ γ br , yet the synchrotron spectrum, ϵL ϵ , peaks at an energy that is defined by the proton Lorentz factor.Depending on the slope of the target-photon population, the synchrotron energy spectrum is dominated by the emission of the highest energy protons, i.e. ∝ B ′ γ 2 p,max (e.g. for s ph = 2 most target-photons are of low energy providing threshold targets for the highest energy protons of the distribution), or the emission of the lowest energy protons, i.e. ∝ B ′ γ 2 p,min (e.g. for s ph = −1 most of the target-photons have high energies providing threshold targets for the low-energy protons).

General remarks
Before closing this section we summarize our main findings about the Bethe-Heitler pair distributions.
For monoenergetic protons interacting with monoenergetic photons we find that: • most pairs are produced with a Lorentz factor ∼ ϵ −1 ; this was also pointed out by [26].
• the injection rate reaches a peak value when the protons interact with target photons near the threshold.Thus, the pair population will be dominated by near-threshold interactions.
• most of the injected energy is transferred to pairs that have a Lorentz factor equal to that of the proton, i.e. γ e ≈ γ p , for near-threshold interactions.In far-from-threshold interactions roughly equal amounts of energy per logarithmic decade are transferred to pairs.
For PL protons distributions interacting with a PL photon distribution, we find that (see figure 3): • the produced pair energy spectrum, γ 2 e Q(γ), is determined by the interactions of the most energetic protons of the distribution, when these dominate in number (i.e., s p ≤ 1).More specifically, Figure 4: Steady-state number distribution of Bethe-Heitler pairs and their emitted synchrotron spectrum for the same parameters as those used in figure 2 and B ′ = 40 G.Each curve, starting from the bluer one, is multiplied by a factor of [1,4,16,1,4] respectively to avoid overlap.While the cooling Lorentz factor determines the peak of the steady-state pair distribution regardless of s ph , the peak energy of the emitted synchrotron energy spectra depends on the target photon index.
when most of the PL photons are interacting far from threshold with the high-energy protons, then the pair energy distribution is determined by the away-from threshold interactions resulting in an almost flat distribution.when most of the PL photons are interacting near the threshold with the high-energy protons, then the shape of the pair energy distribution is determined by interactions close to the threshold, and is described by a PL component with a well defined peak at γ e ≈ γ p,max .
• when the lowest energy protons of the distribution dominate in number (s p ≥ 1), the produced pair energy spectrum is: determined by the low-energy proton produced pairs when s p − s ph ≳ 1, determined by the high-energy proton produced pairs when s p − s ph ≲ 1, almost flat when s p − s ph ≃ 1.
In astrophysical sources, such as AGN, the target photon fields are typically soft distributions with s ph ≳ 1, originating by non-thermal mechanisms that cover a wide energy range starting e.g. from a few eV and extending to a few MeV.These extended photon fields provide near-threshold targets for protons of different Lorentz factors.As a result the electromagnetic component originating by Bethe-Heitler synchrotron emission is expected to have a peak value determined by those protons of the distribution that carry most of the energy.

Bethe-Heitler γ-ray emission in blazars
Various processes have been invoked to explain the γ-ray emission from blazars, including electron inverse Compton scattering on synchrotron photons [e.g.[29][30][31] Table 3: Typical values of the peak frequencies (and corresponding photon energies in units of m e c 2 ) of the two SED humps for LSPs, ISPs and HSPs [38].

Analytical considerations
We present analytical approximations for the locations and the luminosities of the two spectral peaks as they appear in a plot of log 10 (ν F(ν)) versus log 10 (ν).We also discuss the relative importance of photopion interactions and photon-photon pair production.We assume that all particle populations are effectively monoenergetic, namely they have distributions that are energetically dominated by particles of a characteristic energy.Moreover, a particle with a given Lorentz factor is assumed to produce photons at the characteristic synchrotron frequency [28].

Characteristic photon energies
The low-energy peak of the SED appears at an observed frequency ν S (see Table 3), which is related to the characteristic Lorentz factor, γ ′ e , of accelerated electrons in the jet (henceforth, primaries) as is a dimensionless measure of the co-moving magnetic field strength B ′ .Similarly, the high-energy peak of the SED, which appears at a frequency ν γ (see Table 3), is created by the synchrotron radiation of Bethe-Heitler pairs with characteristic Lorentz factor γ ′ BH , In section 2 we showed that the energy distribution of pairs created via Bethe-Heitler pair production peaks at a Lorentz factor almost equal to the Lorentz factor of the parent proton, i.e. γ ′ BH ≈ γ ′ p .These are the pairs that will contribute to the peak of the emitted synchrotron spectrum (see section 2.2).Therefore, the peak of the high-energy component will appear at The ratio of the SED peak energies is then simply given by .Protons interact at the threshold with electron synchrotron photons, and produce pairs that emit synchrotron radiation within the observed γ-ray limits (see Table 3) for γ ′ e and γ ′ p values selected from the intersection of colored lines with the shaded region for each blazar subclass.The peak energy of the SSC component is lower than the peak energy of the Bethe-Heitler synchrotron component for parameter values drawn from the magenta shaded region.There are combinations of γ ′ p and γ ′ e leading to Bethe-Heitler synchrotron spectra peaking in γ rays.
similarly to the proton-synchrotron scenario but without the multiplying factor m e /m p [35,39].Moreover, the primary electron and proton Lorentz factors are not always independent model parameters, like in the proton synchrotron model.More specifically, assuming that protons with Lorentz factor γ ′ p interact with synchrotron photons of energy ϵ ′ S at the threshold for Bethe-Heitler pair production, then the following condition has to be met, Using Eqs.(3.4) and (3.5) we create a parameter space for γ ′ e and γ ′ p for each source class and three indicative magnetic field strengths -see figure 5. Colored lines indicate the locus of points satisfying the energy threshold condition of Eq. (3.5), while the shaded-grey region corresponds to Eq. (3.4) for a range of ϵ S and ϵ γ values that are representative for each blazar class.For a specific magnetic field strength and source class, there is small range of γ ′ e and γ ′ p values that satisfy both conditions simultaneously.
Moreover, by combining Eqs.(3.3) and (3.5) we can express the dimensionless energy of the high-energy hump as where we introduced the notation f X ≡ f /10 X (in cgs units, unless stated otherwise).Interestingly, the peak energy of the proton synchrotron spectrum in this scenario is a constant fraction of ϵ γ , namely The synchrotron self-Compton (SSC) emission produced by primary electrons peaks at a photon energy where we assumed that the scatterings take place in the Thomson regime.The SSC component peaks at lower energies than the proton synchrotron spectrum, see Eq. (3.7), when Additionally, the SSC component peaks below the Bethe-Heitler spectrum, if This condition is satisfied for most parameters that lead to a Bethe-Heitler component peaking in γ-rays, as shown in figure 5.
An isotropic external thermal photon field, if present, can also provide Bethe-Heitler threshold targets for relativistic protons in the jet.We generally describe the thermal emission with a grey body of temperature T in the AGN rest frame (e.g.infrared emission from dusty torus).We also assume that the jet emitting region (blob) is moving within the external radiation field, which is considered isotropic in the AGN rest frame.Hence, the temperature and energy density of the external photon field, as measured in the blob rest frame, are relativistically boosted as T ′ = ΓT and u ′ ext = Γ 2 u ext [40].The threshold condition then reads where we considered a typical synchrotron energy for LSP blazars.

Peak photon luminosities
We continue with an estimation of the peak luminosities of the SED humps.Under the assumption of monoenergetic particle population and the δ-function approximation for the synchrotron emissivity, the low-energy peak luminosity can be written as where u ′ B = B ′2 /(8π), N e is the number of primary electrons with Lorentz factor γ ′ e in the blob, and the factor in square brackets accounts for the cooling of monoenergetic electrons.The number N e can be related to the electron energy density u ′ e as At this point it is useful to introduce the electron compactness which is a dimensionless measure of the total energy density and an input parameter of the numerical code ATHEνA, results of which will be presented in section 3.3.By combining Eqs.(3.13), (3.14) and (3.15), we obtain the low-energy peak luminosity (in the observer's frame) In order to compute the peak synchrotron luminosity emitted by Bethe-Heitler pairs, we will assume that the latter radiate all of their energy via synchrotron, namely the injected power in Bethe-Heitler pairs equals the power emitted in synchrotron photons L ′ BH,S = L ′ BH,e .The total energy transferred to Bethe-Heitler pairs per unit time is a fraction of the relativistic proton luminosity L ′ p , where f BH ≤ 1 is the Bethe-Heitler pair production efficiency, σBH ≃ 8 × 10 −31 cm 2 is the maximum effective cross section accounting for the inelasticity of the interaction [41], and n ′ t is the target photon number density.If ϵ ′ t is the characteristic target photon energy, then the number density can be written as where in the last expression we introduced the photon compactness ℓ t , defined in a similar manner as in Eq. (3.15).We can also express L ′ p in terms of the proton compactness ℓ p as Combining the equation above with Eqs.(3.17) and (3.18) we obtain the final expression for the peak luminosity of the high-energy component of the SED (in the observer's frame) In the special case where the target photons for Bethe-Heitler pair production are the primary electron synchrotron photons, which may be relevant for HSP blazars, then ϵ ′ t = ϵ S (1 + z)/δ and the photon compactness is given by where we also used Eq.(3.13).The ratio of the peak luminosities, which is the equivalent of the Compton ratio in purely leptonic models, can then be written as If the target photons were provided by an external source of radiation of energy density u ′ t ≃ Γ 2 u ext and photon energy ϵ ′ t = Γϵ ext , then the ratio of the peak luminosities would read (3.23) In the proposed scenario the peak luminosity ratio has a different dependence on ϵ S depending on the dominant target photon field for Bethe-Heitler pair production: synchrotron jet emission or external radiation.By varying model parameters, like ℓ p and ℓ e , it is possible to obtain a wide range of values for the peak luminosity ratio.We will compare the predictions of the fully numerical models against the distribution of observed values in Sec.3.3.

Photopion and photon-photon pair production processes
Relativistic protons may pion produce on low-energy photons, leading eventually to the production of energetic γ-ray photons, pairs, and neutrinos.It is therefore interesting to compare the expected luminosities of synchrotron-emitting pairs and neutrinos with the γ-ray luminosity of Bethe-Heitler pairs.Moreover, we discuss the role of photon-photon (γγ) pair production, which is relevant for the attenuation of energetic photons from pion decay.
For simplicity, we assume that photopion interactions take place at the ∆ + resonance, where π 0 and π + pions are produced at a ratio 1 : 1 (i.e.equally probable).Neutral pions will decay into two photons, while charged pions will produce positrons/electrons and neutrinos, after a series of reactions.We consider that the energy lost by the reactants is equally distributed to the products, resulting in γ-ray photons from π 0 decay with energy, where we used that the inelasticity of p-γ interactions is κ pγ ≃ 0.2 close to the ∆ + resonance.Similarly, the neutrinos produced in the π + decay chain will acquire 3/4 of the pion energy, leaving 1/4 of the latter for the positron created.The energy injected per unit time in γ-ray photons and neutrinos of all flavors can be then estimated as: where R ′ b is the blob radius, L ′ p the proton population luminosity, κ pγ σ pγ = σpγ ≃ 7 • 10 −29 cm 2 is the product of the inelasticity and the cross section of the photopion interaction [42] and n ′ t,pγ is the number density of the relevant target photons.
We compare next the luminosities of secondaries from photopion interactions with the luminosity of synchrotron photons emitted by Bethe-Heitler pairs, see Eq. (3.17) In the case of monoenergetic protons, as assumed here for simplicity, the luminosity ratio depends only on the cross sections and number densities of the target photons of the interaction processes.In AGN jets where photons are produced by non thermal mechanisms, the differential number density of photons can be usually described by a power law n ′ t ∝ ϵ ′−s ph +1 .Therefore, the number density of target photons for each interaction will generally depend on the photon index s ph .In the special case of s ph = 1, where the photons are equally distributed in number across the energy spectrum, we expect n ′ t,pγ /n ′ t,BH = 1 and L ν /L BH,S ≈ 30.For s ph 1 the target number ratio would be n ′ t,pγ /n ′ t,BH ≈ 10 2(−s ph +1) , since ϵ pγ,thr /ϵ BH,thr ≈ 140.For example, if s ph = 2, then n ′ t,pγ /n ′ t,BH ≈ 0.01 resulting in L ν /L BH,S ≈ 0.1.Summarizing, for photon indices in the range 1 to 2, we could expect 0.1 ≲ L ν /L BH,S ≲ 10.
The energetic photons produced by the decay of neutral pions are expected to be attenuated by low-energy photons in the source.We can estimate the opacity of the source to γγ pair production of energetic photons with observed energy as ϵ γ = δϵ ′ γ /(1 + z) as where we used the approximate γγ cross section from Ref. [43].In figure 6 we present indicative results for LSP, ISP, and HSP blazars.The optical depth for each case is integrated over a broadband target photon distribution 3 .Figure 6 shows the transition from an optically thin source (τ γγ ≤ 1) to an optically thick source (τ γγ ≥ 1).When the latter transition happens we can observe that the luminosity of the electromagnetic spectrum drops significantly due to the internal photon-photon annihilation that becomes important.For the examples considered, the source is optically thick to energetic γ-rays from pion decay.

Numerical approach
According to the analytical results presented in the previous section, it is possible to find parameter values that can lead to a high-energy SED component powered by synchrotron radiation of Bethe-Heitler pairs.This analysis is however simplified, because it considers effectively monoenergetic particle distributions.Moreover, the effects of other processes, like photomeson interactions and photon-photon pair production were evaluated in an approximate way.
In order to examine in more detail the hypothesis that Bethe-Heitler pair production can produce the γ-ray component in blazar spectra, we will use the leptohadronic code ATHEνA [18,24].The code numerically solves the kinetic equations that describe the evolution of stable particle distributions that are injected or created inside a spherical magnetized source.More specifically, the code solves a system of coupled partial differential equations for protons, neutrons, electron-positron pairs, photons, and neutrinos.Unstable particles, such as kaons, pions, and muons, are assumed to decay instantly into secondary lighter particles, hence only their injection rates are computed.All stable particles are also allowed to escape the source on a timescale equal to the light crossing time R ′ b /c.The physical processes that are included are: synchrotron emission from charged particles, synchrotron self-absorption, electron inverse Compton scattering, photon-photon pair production, and proton-photon interactions that may lead to pion production (photomeson) or pair production (Bethe-Heitler) -for the numerical implementation, see [18].
The code takes as input the parameters describing the emitting region, namely radius and magnetic field strength, and the characteristics of the power-law distributions of electrons and protons injected in the blob (i.e.minimum and maximum Lorentz factors, slope and compactness).Given these input parameters, we evolve the system long enough to reach a steady state (i.e.∼ 10 R ′ b /c for most cases).Based on the steady-state photon number density, we compute the escaping photon flux in the observer's frame after performing the appropriate Doppler boosting [44].We numerically investigate the parameter space of the model in search for parameters that produce a high-energy SED component dominated by synchrotron emission of Bethe-Heitler pairs for different blazar classes (LSPs, ISPS, and HSPs).As indicative values for the peak photon energy of each subclass we use those listed in Table 3.
We start by selecting 3 indicative values for the blob radius, R ′ b ∈ {10 14 , 10 15 , 10 16 } cm.For each value we then use 3 magnetic field strengths, B ′ ∈ {0.5, 5, 50, 500} G, in order to cover a wide range of physically plausible values.Then, for each pair of (R ′ b , B ′ ) values and for each spectral subclass with fixed (ϵ S , ϵ γ ), we determine γ ′ e and γ ′ p using the analytical expressions of the previous section as a guide (see e.g.Eqs.3.3, 3.5 and 3.11 and figure 5).More precisely, the distributions of injected particles are modeled as power laws that energetically peak close to the analytically predicted values.We then vary the normalization of the particle distributions, or equivalently ℓ e and ℓ p , to produce a range of peak luminosities that fall within the range of observed values.
For LSPs we also consider the presence of an external thermal photon field with a grey-body spectrum.The temperature is chosen so that the thermal photons from the peak of the distribution fulfil the Bethe-Heitler threshold condition in the blob rest frame, see Eq. (3.11).The compactness (or comoving density) of the external photon field is adjusted so that it leads to values of A (ext)  γS > 1, as usually observed in LSP spectra [38].The plausibility of the chosen values will be discussed in section 3.3.3.In figure 7 we present a ISP-like blazar spectrum with different spectral components highlighted (the parameters used are the same as those listed in Table 4 for the ISP case with B ′ = 50 G and R ′ b = 10 14 cm).Starting from low energies, we identify first the synchrotron spectrum of primary electrons, Their distribution is steeper than the injected power law due to synchrotron cooling.The analytical expectation for the peak of the synchrotron spectrum in this case would be ∝ Bγ 2 e,max , but this overestimates the position of the numerical spectrum by a factor of ∼ 7. The reason for the discrepancy is that the cooled electron distribution falls off faster than the asymptotic power law with −s e − 1 due to the narrow dynamic range of the injected power law.As a result, the electrons that contribute most to the peak of the radiated synchrotron spectrum have γ ∼ γ e,max /2.5.The primary synchrotron component is followed by a subdominant, in this example, synchrotron self-Compton bump (dash-dotted blue line); its peak energy agrees with the estimate of Eq. (3.8).At slightly higher frequencies, we observe a distinct bump due to proton synchrotron radiation (dashed red curve), which peaks at energy given by Eq. (3.7).At ∼ 10 − 100 GeV energies we observe the synchrotron emission of Bethe-Heitler produced pairs (thick dashed yellow line) when γγ absorption is neglected; the spectrum is computed using the analytical function for the pair injection, see Eq. (2.2).A two-component spectrum (dotted magenta line) emerges at the highest energies, which consists of synchrotron emission of pairs from charged pion decays and photon emission from neu-tral pion decay.These energetic photons are attenuated, however, inside the source (see e.g.yellow lines in figure 6), producing electron-positron pairs that, for the specific parameters, emit synchrotron photons across a wide range of energies between 10 keV and 100 GeV.The synchrotron spectrum of the γγ-produced pairs is shown as a separate component marked with cyan color (with attenuation included).For the selected parameters is less luminous than the Bethe-Heitler and proton synchrotron components, but emerges as a separate component between 10 keV and 100 keV.Finally, the black line shows the spectrum after accounting for attenuation due to the extragalactic background light (EBL).For this example, we considered a source at redshift z = 0.3 and used the EBL model of Franceschini [45].Due to the combined attenuation, only part of the Bethe-Heitler component would be visible.The magenta line in the high-energy part of the spectrum describes the energy distribution of neutrinos and anti-neutrinos of all flavors.The peak neutrino luminosity is comparable to the peak luminosity of the Bethe-Heitler component before attenuation, a result that agrees with Eq. (3.26).

Blazar SEDs
In figure 8 we present an ensemble of theoretical SEDs that were built as outlined in Sec.3.2, assuming an emitting region with radius R ′ b = 10 15 cm; all other parameter values are listed in Table 4 of appendix D. We find parameter sets for all blazar subclasses resulting in a dominant BH-synchrotron component in the γ-ray regime (see thick solid lines in each panel) while producing broadband spectra with similar features as the observed ones (e.g.peak photon energies and luminosities).Changes in ℓ e and/or ℓ p are commonly invoked to reproduce flux variability [e.g.[46][47][48].Therefore, for each case, we demonstrate the effect that different values of ℓ e (dotted lines) and ℓ p (dash-dotted lines) have on the SED morphology.In all panels, we also include the expected all-flavor neutrino spectrum for completeness.
In the first row of figure 8 we show spectra that resemble those of LSPs.The high-energy peak remains unaffected by changes in ℓ e because protons pair produce close to the threshold with external photons that have higher number densities than the electron synchrotron photons in the comoving frame.Changes in ℓ p , however, have a linear effect on the peak luminosity of the Bethe-Heitler synchrotron spectrum as long as the target density remains fixed (see Eq. 3.20).As ℓ p increases, the third bump peaking at ∼ 0.1 TeV, which originates from pions created by the proton population interacting with the low-energy photons of the electromagnetic spectrum, becomes progressively more luminous that the Bethe-Heitler component (see cases for B ′ = 5 G and B ′ = 50 G).However, in smaller sources (e.g.R ′ b = 10 14 cm) with stronger magnetic fields, such as B ′ = 50 − 500 G, the same increase in ℓ p leads to drastically different spectra that are almost featureless, characterized by a MeV peak with much higher luminosity (see e.g.dash-dotted lines in middle and right panels of the top row in figure 23).These spectra are characteristic of a non-linear cascade developed in the source as protons lose efficiently their energy through photo-hadronic interactions on target photons that are not any more fixed, but are indirectly related to the proton population: attenuation of energetic photons produced in photo-hadronic interactions leads to the injection of relativistic pairs that produce low-energy target photons [for more details on hadronic supercriticalities see 49,50].Larger sources, like those presented in figures 8 and 24 are less compact for the adopted values, hence we do not observe featureless spectra peaking in MeV energies (with the exception of the LSP-like case with B ′ = 50 G in figure 8).Lower ℓ p values still produce a two-component spectrum, but with lower γ-ray luminosity than the default spectra displayed with solid curves.In this case, higher bolometric luminosities may be achieved for higher values of the Doppler factor.
In the ISP-and HSP-like cases, the primary electron synchrotron photons are the main targets for Bethe-Heitler pair production, thus changing the dependence of the BH-synchrotron component on ℓ e .The target photon density is proportional to ℓ e while the electron inverse Compton luminosity  4).In all panels, dashed lines represent the synchrotron spectrum of Bethe-Heitler pairs computed using the analytical function for the pair injection, see Eq. (2.2).Dash-dotted and dotted lines represent spectra obtained by varying the proton and electron compactness, respectively, by two orders of magnitude with respect to their nominal values (solid lines).The energies of the low-and high-energy humps for the baseline model (solid lines) are indicated with symbols.Long-dashed lines in the LSP graphs represent the external photon field.Gray shaded regions indicate the observational range of peak frequencies for the low-and high-energy SED components.Black horizontal lines at ∼ 1 − 10 PeV represent the neutrino flux prediction according to Eq. (3.26).scales as ℓ 2 e .Therefore, there are cases where an increase of ℓ e can produce an inverse Compton component that is more luminous than the Bethe-Heitler synchrotron -see e.g.upper dotted curves for B ′ = 0.5 G and B ′ = 5 G.We also note that the peak energy of the SSC component is typically lower than the Bethe-Heitler synchrotron peak energy (see Eqs. (3.7) and (3.8)), which was chosen to fall within the typical range of ϵ γ values for ISP and HSPs.Lower values of ℓ e value will result in the suppression of inverse Compton emission from primary electrons.Still, a lower limit on ℓ e is needed to produce a Bethe-Heitler synchrotron component that is more luminous than the proton synchrotron, contrary to the cases displayed by the lower dotted curves in the panels for ISPs and HSPs with B ′ = 50 G and B ′ = 500 G.For the ISP and HSP cases we examined, we find a linear dependence of the Bethe-Heitler synchrotron emission on ℓ p .Because of the lower values of ℓ p used to produce the ISP/HSP-like spectra, the system does not enter into the non-linear regime by increasing ℓ p , as shown in upper panels for the LSPs of figure 23 in Appendix C. Finally, we note that photon-photon pair production severely attenuates the Bethe-Heitler synchrotron component in most HSP-like cases.As a result, the high-energy SED component does not peak within the range of observed peak energies, especially for cases with B ′ ≥ 5 G.In larger sources (e.g.R ′ b = 10 15 cm and R ′ b = 10 16 cm) with B ′ = 0.5 G the photon density of the synchrotron-produced photons is lower and the attenuation due to γγ pair production not being as intense as for stronger magnetic fields.
In figure 8 we also observe that the neutrino luminosity is in most cases comparable to the Bethe-Heitler synchrotron peak luminosity, since the targets for BH pair production are comparable to those for photopion production.In such cases the analytical prediction of Eq. (3.26) states that the Bethe-Heitler synchrotron and the photopion-related luminosities will differ from each other no more than an order of magnitude.In ISP-and HSP-like models we require the primary electron synchrotron photons to provide on-threshold targets for Bethe-Heitler interactions of protons.However, the on-threshold target photons for photomeson interactions are ∼ 150 times more energetic than the threshold photons for Bethe-Heitler interactions.In some models the photomeson and Bethe-Heitler target photons are both part of the power-law synchrotron spectrum of primary electrons.In these cases (e.g.HSP-like models in figure 8), the luminosity ratio of Bethe-Heitler synchrotron photons and neutrinos may range from 0.1 to 10 depending on the photon index of the target photon spectrum (see also section 3.1.3).Otherwise, if the photomeson targets are part of the exponential cut-off of the electron synchrotron spectrum, their number density is much smaller than the density of photons responsible for Bethe-Heitler pair production, resulting in a much lower neutrino luminosity than the one of the Bethe-Heitler component (e.g.ISP-like cases in figure 8).
At this point it would be interesting to compare the accuracy of the analytical estimations presented in section 3.1 about the injected and synchrotron emitted luminosities of primary electrons and Bethe-Heitler pairs against our numerical results that are not limited by simplifying assumptions.Figure 9 shows the ratio of the synchrotron radiated luminosity over the injected one for both the primary and secondary electrons (panel a), the ratio of the numerically calculated synchrotron luminosity over the analytical expectation for the primary electrons (panel b), and the same ratio for the Bethe-Heitler pairs (panel c).In all panels, the ratios are plotted against t syn /t dyn , where t dyn = R ′ b /c, which can be considered as a proxy for the synchrotron cooling efficiency of particles injected with Lorentz factor γ e .Panel (a) shows that Bethe-Heitler pairs are fast cooling due to synchrotron radiation for all models, since t syn /t dyn ≤ 1, while for the primary electrons this is true only for models with B ′ ≥ 50 G (regardless of the source size).Our analytical formulas (see e.g., Eqs.3.16 and 3.20) were derived assuming that all the injected luminosity of the electrons is radiated through synchrotron, which is not always the case for the primary population.As a result, we expect that our analytical predictions for the primary electrons will not match the numerical results when electrons do not cool efficiently due to synchrotron (i.e., t syn /t dyn > 1).The latter can, indeed, be observed in panel (b) where the ratio of the numerical and the predicted synchrotron luminosities is less than unity for models with B ′ ≤ 5 G (i.e., the analytical expression in this regime overestimates the primary electron synchrotron luminosity).
To analytically calculate the Bethe-Heitler synchrotron luminosity (Eq.3.20) we assumed that (i) pairs radiate all their energy due to synchrotron and (ii) the target photons for pair production in ISP-and HSP-like case are the primary electron synchrotron photons with a luminosity given by 3.16.While the first assumption is always fulfilled, as shown in panel (a) of figure 9, the second is not always true; the primary synchrotron luminosity is not always equal to the theoretical value given by 3.16 due to inefficient synchrotron cooling, as shown by panel (b) (black and blue markers).For the LSP-like cases the dominant targets for Bethe-Heitler pair production are assumed to be of external origin with a thermal energy distribution.In these cases, the available targets might also be of lower luminosity value than the theoretical one, being the grey body bolometric luminosity.However, since the protons interact with the near-peak photons of the spectrum and not with the low-energy tail of the grey body distribution, not all the grey body luminosity is available for pair creation.Panel (c) shows that Eq. 3.20 is good proxy of the Bethe-Heitler synchrotron luminosity (within factors of 2-3) when the primary electron luminosity is radiated through synchrotron (red and yellow markers), but it can overestimate by a factor of 10-30 the luminosity when primary electrons are not fast synchrotron cooling (blue and black markers).

Compton Dominance
Next we investigate how our theoretical SEDs compare to observations of Fermi-detected blazars.
In addition to the SED models shown in figure 8 we compute photon spectra for R ′ b = 10 14 cm and 10 16 cm (see figures 23 and 24, and Table 4 for a complete list of parameter values).For each model we calculate the peak synchrotron luminosity and peak frequency (L pk,S and ν S ), and the peak luminosity L pk,γ of the high-energy component.We then define the Compton dominance as the ratio of the two peak luminosities, A γS = L pk,γ /L pk,S .In figure 10 we compare our results against those derived from a sample of 781 Fermi-detected AGN, which includes 504 FSRQs and 277 BL Lac objects [51].[51].The size of colored markers indicates the size of the blob, with the smallest markers corresponding to 10 14 cm and the largest to 10 16 cm.Open markers indicate models that (i) fall outside the locus of observed points, or (ii) have γ-ray peak frequencies outside the typical range of values for each subclass (see table 3), or (iii) require non-physical combinations of temperature and luminosity of the external photon field (see section 3.3.3).
Most theoretical models fall within the range of observational values and follow the same trend as the latter.The scatter in our models is attributed to different sizes of the emitting region and Doppler factor values.While this agreement with the data is achieved partially by construction (i.e., we selected our parameter values in order to reproduce peak luminosities and frequencies, as explained in section 3), it was not clear a priori if the Bethe-Heitler emission from pairs could dominate the high-energy emission in different blazar subclasses.This might be an indication that systems where Bethe-Heitler pairs play an important part in the production of the high-energy component of the blazar SED can exist and be hidden amongst the already observed AGN spectra.
However, some models can be rejected if they do not pass certain criteria, which we describe below.Criterion (i): if a model falls well outside the observational locus of points in the peak γ-ray luminosity (or Compton dominance) versus peak synchrotron frequency diagrams, then we characterize it as non-plausible (e.g., ISP-like model with B ′ = 500 G and R ′ b = 10 15 cm).Criterion (ii): models might not qualify as plausible candidates for a certain blazar subclass, if the peak frequency of the high-energy hump does not fall into the typical range of values for that subclass (see table 3 and grey-shaded areas in figure 8).Criterion (iii): LSP-like models that require non-physical combinations of temperature and luminosity for the external photon field are also excluded (for more details see the next section).Rejected models are indicated with open markers in figure 10.
Figure 11 shows both the total number of models we computed, 12 for each type (LSP-, ISP-, HSP-like), and the number of the accepted models based on the criteria described in the previous paragraph.In summary, we found that almost all HSP-like models (9/12) can be excluded solely by criterion (ii).Most LSP-like models (7/12) can be excluded based on criterion (iii), while ISP-like models (2/12) can be excluded either on criterion (i) or (ii).Our results suggest that if we were to observe blazars with a Bethe-Heitler-dominated γ-ray spectrum, most of them (55%) would be ISPs, a few LSPs (28%) and, only a a minority would be HSPs (17%).

External photon fields
To reproduce LSP-like spectra within our scenario we invoked arbitrary external radiation fields to compensate for the low energy of jet synchrotron photons.More specifically, these photons are not energetic enough (for typical values in LSPs) to provide near-threshold targets for Bethe-Heitler interactions, which are needed to produce a luminous high-energy component.Therefore, a greybody photon field with effective temperature T ′ ≈ 10 2 − 10 5 K in the blob frame was introduced.Considering a Doppler factor of δ ∈ [40 − 80] and Γ = δ/2 (see table 4 in appendix D), the aforementioned temperatures translate to T = T ′ /Γ ≈ 10 − 10 3 K in the AGN rest frame.The energy density of the thermal field (in the AGN rest frame) is related to the photon compactness ℓ ext as u ext = 3ℓ ext m e c 2 /Γ 2 σ T R ′ b .By requiring u ext < u BB , where u BB = αT 4 is the energy density of a black body with effective temperature T and α is the radiation density constant, we can set an upper limit on ℓ ext , In most LSP-like sources, the invoked thermal field does not satisfy the limit above (see table 4 in appendix D), making this a non physical choice.An attempt to lower ℓ ext , while trying to produce the same L γ through Bethe-Heitler synchrotron radiation, would require a higher proton compactness (see Eq. 3.20).This, in turn, would increase the energetic requirements of the model (see also section 3.3.4)and change the relative ratio of the proton synchrotron and pair synchrotron peak luminosities, leading to dual-component γ-ray spectra (see e.g.solid red line in the top right panel of figure 8).From this analysis, only LSP-like models with R ′ = 10 15 − 10 16 cm and B ′ = 5 − 500 G are viable.

Jet power
The power of a two-sided blazar jet, ignoring the contribution of cold protons, can be expressed as [e.g.[52][53][54]: i=B,e,p,ph where R ′ b is the radius of the region that produces the steady emission, and is assumed to be equal to the cross-sectional radius of the jet, Γ is the jet Lorentz factor, β ≈ 1 is the jet speed (in units of c), u ′ i and P ′ i = 1 3 u ′ i are the energy density and the pressure of the magnetic field and of the relativistic particles in the emitting region.In figure 12 we plot the jet power of our models against the proton-to-magnetic energy density ratio (left panel) and the γ-ray peak luminosity (right panel).Note that the jet Lorentz factor is not the same for all cases (see Table 4).Models with large emitting regions and strong magnetic fields are characterized by high powers that surpass the Eddington luminosity for a 10 8 M ⊙ black hole mass (horizontal line) by several orders of magnitude.While the Eddington luminosity of the accreting black hole is not a strict physical limit to the power of the jet, such high values are not expected in sub-Eddington accreting black holes.The jet power in an accreting system can be written as P jet = η j Ṁc 2 , where Ṁ is the accretion rate onto the black hole and η j is the jet-formation efficiency, which can be ∼ 1.5 at most for magnetically arrested accretion disks [55,56].Thus, the jet power can be as high as P jet ≲ 15 ṀEdd c 2 /η −1 ≃ 1.5 × 10 47 M 8 /η −1 erg s −1 , where η is the radiative efficiency of the accretion flow.Furthermore, the displayed powers are actually lower limits (due to the assumed very narrow proton energy distribution and the neglect of the non-relativistic plasma component).Therefore, models with P jet ≫ 10 47 erg s −1 can be discarded based on energetic grounds.Moreover, blazar jets are thought to be strongly magnetized and narrow near their base, with weaker magnetic fields on larger scales [57,58].It is therefore difficult to physically explain the formation of large emitting regions (at pc scales) with hundreds of Gauss (comoving) magnetic field strengths.These arguments also disfavor the models with extreme jet powers.From this analysis, the physically plausible models are those with: R ′ b = 10 14 cm and B ′ ≤ 100 G, R ′ b = 10 15 cm and B ′ < 100 G, and R ′ b = 10 16 cm with B ′ ≲ 5 G.

Summary and Discussion
In this paper we took a closer look at the characteristics of the distribution of relativistic electrons and positrons produced in Bethe-Heitler interactions of relativistic protons with low-energy photons.
While most pairs are created with a Lorentz factor γ e ≈ ϵ −1 , most of the energy is transferred to particles with γ e ≈ γ p for near-threshold interactions between protons and photons of single energy ϵ.Far-from-threshold interactions lead to more extended pair distributions, and tend to distribute almost equal amounts of energy per logarithmic decade in γ e .In astrophysical environments where interactions of protons take typically place on extended target photon distributions, we find that the energy distribution of pairs is usually determined by near-threshold interactions, and has a welldefined peak energy.Additionally, we provide an empirical function that is able to describe the Bethe-Heitler pair distribution satisfactorily.The empirical function can be implemented in numerical codes and replace the time-consuming calculation of double integrals.Our function, however, has not been constructed to model pair distributions for interaction energies γ p ϵ < 2. The contribution of these interactions to the total pair production spectrum has been found to be smaller by an order of magnitude than the contribution of interactions with γ p ϵ ≈ 5 (very close to the threshold, the effective cross section of the process is very small -see e.g.figure 1 in [26]).As a result, when target photons for Bethe-Heitler are described by an extended distribution, γ p ϵ < 2 interactions will not be important in shaping the overall secondary pair distribution.
Furthermore, we examined if it is possible to produce γ-rays from energetic Bethe-Heitler pairs in blazar jets, and produce SEDs that resemble the observed ones.We used analytical arguments (under several simplifying assumptions) to determine the relevant parameter space, and created 12 numerical SED models for each blazar subclass (36 in total) using the code ATHEνA.As the numerical models do not suffer from the simplifying assumptions used in our analytical estimates, they were also used to cross check the validity of the latter.We found that the parameter space leading to a dominant Bethe-Heitler γ-ray component is very constrained, especially when the jet photons are the targets for Bethe-Heitler interactions (see e.g.figure 5).Moreover, small changes in parameters like the energy injected to accelerated electrons and protons, can favor other processes, like SSC or proton-synchrotron radiation, and eventually hide the Bethe-Heitler γ-ray emission (figure 23).The latter attribute can motivate small adjustments of the injected proton and electron luminosities in order to fine tune the proton synchrotron and Bethe-Heitler synchrotron peak luminosities.There are cases (see e.g.figure 7) where the aforementioned components are both distinguishable with comparable luminosities.But by changing the injection luminosities by a factor of 0.1-0.2 in logarithm, one could suppress the proton synchrotron emission and acquire a spectrum with a broad γ-ray component, the Bethe-Heitler synchrotron one.Given the limited number of models we created, we did not use the presence of additional spectral components in the broadband SED to exclude said models.However, we implemented three other criteria to determine whether a model is accepted as a possible candidate for a blazar subclass or not.First, we checked if the peak position of the γ-ray component fell inside the frequency range set in table 3 for each blazar type.Second, we checked whether the numerical models fell within the locus of observed points in figure 10.Lastly, for LSP-like models, we checked if the combination of temperature and density for the external photon field was physical or violated the black body limit.Interestingly, after enforcing the three criteria, we found that most HSP-like models can be excluded because of the strong attenuation of γ-rays above a few tens of GeV (see e.g.figures 6 and 23).Nevertheless, a few LSP-like and most ISP-like models passed all criteria, suggesting that sources with Bethe-Heitler γ-ray emission could exist in the blazar population, but they would not be the most common ones.Interestingly, in an independent work where SED fitting was performed for a sample of 34 blazars, the best-fit model for 3 ISPs yielded a dominant Bethe-Heitler synchrotron component in γ-rays (Rodrigues, X. et al., in preparation).
Our blazar models were computed using narrow power-law proton distributions, even though acceleration mechanisms create wide distributions, typically starting from the proton rest mass energy.Therefore, questions about the contribution of lower energy protons to the total spectra and jet power may arise.If the proton population is hard (s p < 2), we find that the low-energy protons do not significantly contribute to Bethe-Heitler pair production and the total jet power, since most of the proton population energy is carried by protons at the high-energy tail of the distribution.These energetic protons are responsible for the production of pairs emitting in γ-rays.When we extend the proton distributions down to γ ′ p ≈ 1, the compactness is slightly affected (by ∼ 0.7 in logarithm).For soft power laws (s p > 2) extending down to γ ′ p ∼ 1, the required jet power would be order of magnitude higher than the values displayed in figure 12, to compensate for the "inactive" lower energy protons of the distribution that would not contribute to pair production.Alternatively, protons can accelerate into a broken power law with a low energy branch having a slope s p,lo ≤ 2 from γ ′ p ∼ 1 up to γ ′ p,min , and s p > 2 for γ ′ p > γ ′ p,min .When we extend the proton distributions down to γ ′ p ≈ 1, assuming s p,lo = 2 or 1, the compactness is increased by a factor ≈ 18 or 2, respectively.As a result, the main conclusions presented in section 3 remain valid, with the assumption of a proton distribution following a broken power law distribution with a slope of < 2 for γ ′ p < γ ′ p,min .The flat spectrum radio quasar 3C 279 is one of the brightest γ-ray blazars, also classified as an LSP source, located at redshift z = 0.536 [59].In June 2015, the Large Area Telescope (LAT) on board of the Fermi satellite recorded a luminous γ-ray flare in GeV energies (L γ ∼ 10 49 erg s −1 ) with 5-minute timescale variability [60].The short variability timescale, the large amount of energy released in GeV γ-rays, and the large Compton ratio (about a 100) of the flare make it difficult to explain this event with standard leptonic or lepto-hadronic models [for details, see 60, 61].Ackermann et.al [60] proposed an alternative leptonic scenario in which the γ-ray component of the event is produced by a second, energetic electron population, with no specified origin, emitting through synchrotron.This proposal motivated us to investigate whether Bethe-Heitler pair production was able to produce such a relativistic leptonic population.We found that is possible to explain the SED of the 2015 flare if certain conditions are met.First, the shape of the flaring γ-ray spectrum (see figure 4 in [60]) is not very broad around its peak, thus requiring near-threshold Bethe-Heitler interactions to dominate the pair production.Second, a strong magnetic field, e.g.B ′ = 150 G, is needed for boosting the luminosity of the pair synchrotron spectrum.However, even such strong magnetic fields were not enough to reproduce the observed luminosity, so a high Doppler factor value, δ ≈ 190, was needed, as well.For an adopted B ′ value, the proton Lorentz factor can be determined so that the pair distribution peaks at the observed γ-ray energy (see Eq. 3.3).We can then determine the target photon energy that satisfies the Bethe-Heitler threshold condition (see Eq. 3.11).Because of the very high δ values needed, the target photon energy translates to a very low grey body temperature in the AGN rest frame, T ≈ 12 K.Typical temperatures of AGN dust torii are ∼ 300 K to ∼ 1500 K [62].Evidence for colder dust with T∼ 25 K has been reported in blazar Ap Librae, but on hundreds of parsec scales [63].Moreover, the invoked target photon density exceeds that of a black body field with such low temperature.In conclusion, explaining the SED of the 3C 279 2015 flare using a second electron population originating by Bethe-Heitler pair production would require extreme physical conditions, and is therefore deemed unrealistic.
In conclusion, Bethe-Heitler pair production in blazar jets can generate a secondary relativistic lepton population that can emit in γ-rays via synchrotron radiation.Because of the extended proton and photon distributions involved, Bethe-Heitler interactions can happen both near and away from threshold.Nonetheless, their synchrotron spectrum is determined by pairs created by near-threshold interactions.Bethe-Heitler synchrotron radiation can be the dominant γ-ray emission mechanism in low-and intermediate-peaked blazars, but for a limited part of the parameter space.
Each of the three exponential terms appearing in Eq. 2.2 has a specific role in shaping the overall spectrum.First of all, the shape of the energy injection spectrum, namely γ 2 e q BH (γ e ), resembles a logarithmic Gaussian function for interactions near the interaction threshold (see figure 1b), which translates to p(γ p ϵ) ≈ 2 (see figure 1c).We empirically determined that the Gaussian-like shape depends on the interaction energy γ p ϵ.In particular, as the interaction energy increases, the energy injection rate has no longer a well-defined peak but instead becomes flat while extending to lower energies (see e.g.figure 1b).This behavior can be captured by a smaller value of the power-law index p (see figure 14).Since the Gaussian-like function overestimated significantly the energy transferred to the low-energy pairs, a super-exponential term, exp −a 2 2 (γ e,pk /γ e − 1) 2 , was added to make the low-energy part steeper.Finally, the third exponential term controls the high-energy cutoff of the injection spectrum.This, in particular, needed to be softer than the low-energy turnover, hence it is linearly dependent on γ e /γ e,pk .The shape of the three exponential terms of the injection rate function is better illustrated in figure 13.In order to determine the analytical functions for the power-law index p(γ p ϵ) and the normalization A(γ p , ϵ), we used the numerical results presented in figure 1c and 11 additional cases corresponding to different values of the interaction energy γ p ϵ.For each of the 15 benchmark cases we first adjusted the slope p by eye so that the analytical function would match the ATHEνA results.The inferred values of the slope (with a 3% error) are shown with black symbols in figure 14 (left panel).To describe the obtained trend of p with γ p ϵ we introduce the following function, where x ≡ log 10 (γ p ϵ) and a, b, c, s, s 2 , x 0 are free parameters.The first term, ∝ x −s , describes the function near the threshold values of γ p ϵ.It is also multiplied by an exponential term, e −x/x 0 , to suppress its contribution at large interaction energies where the linear term dominates.
To determine the free parameters 4 of the function, we fit the data (black symbols) using emcee [64], a python implementation of the Affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler.This allows a better estimation of the uncertainties and identification of possible degeneracies.We used uniform priors, 100 walkers and propagated each chain for 10,000 steps.Figure 14 (right panel) shows the posterior distributions of the model parameters.The blue dashed curve on the left panel is computed using the median values of model parameters and the gray-shaded area indicates the 68% confidence interval computed from the posterior samples.We also show another solution (solid magenta curve) that is obtained for slightly different parameter values (see black vertical lines in the corner plot).
We next compute the pair injection spectra using the fitted slopes and proceed with a comparison against ATHEνA results -figure 15.The high-energy cutoff of the energy spectrum is extremely sensitive to the value of p as indicated by the large spread in solutions.Therefore, even small changes in the adopted value of p may have a large impact on the energy injection spectra -compare the dashed and solid lines that are obtained using the blue and magenta curves, respectively, shown on the left panel of figure 14.Moreover, the total energy transferred to the pairs better matches the results of ATHEνA when the modified slope fit is used in Eq. (2.2) -see figure 20.For these reasons, we choose to work with the modified slope fit of figure 14 (solid magenta line).After determining the shape of the Bethe-Heitler injection rate function, only the normalization A (see Eq. (2.2)) was left to fit.We first noticed that the logarithm of the peak of the injection spectrum depends linearly on the logarithm of the target photon energy for a fixed proton Lorentz factor (see magenta markers in figure 16).Using the same 15 benchmark cases as in the slope fit, we were also able to determine the dependence of log 10 (A) on the interaction energy, γ p ϵ (see red markers in figure 16).In order to check if the normalization depended also on γ p separately, we used five additional cases that are displayed in figure 1b.These injection spectra where computed for different values of γ p that still corresponded to the same interaction energy γ p ϵ (by appropriately selecting the target photon energy ϵ).The dependence of log 10 (A) on γ p turns out to be the weakest, as shown by the blue markers in figure 14).Based on these trends we devised the function presented in Eq. (2.6), which has 11 free parameters.To determine these parameters we performed a combined fit to the three sets of points displayed in figure 16) using emcee.We employed 500 walkers and each chain was propagated for 5,000 steps.We used uniform prior distributions for all parameters.The corner plot showing the posterior distributions is shown in figure 17, and the grey-shaded bands in figure 16 display the 68% confidence range.Finally, following the same process as we did for the slope p(γ p ϵ), we slightly adjusted the optimal parameters we obtained from the MCMC fitting (see solid black lines in the corner plot) to ensure better agreement of the overall spectral shape with the ATHEνA results .
We finally compare the spectra obtained using the benchmarked/optimized empirical function against the numerical results of ATHEνA.Our results are presented in figures 18 and 19 for the same cases displayed in figure 1.The function given by Eq. (2.2) describes well the peak of the numerical distribution of the Bethe-Heitler produced pairs as a function of the proton Lorentz factor (see left plot of figure 18) and the energy of the target photon (see left plot of figure 19).Similarly, the empirical function describes very well the energy distributions of pairs -see right panels of figures 18 and 19 respectively.

A.2 Integral quantities
In this section we test the performance of the analytical function of Eq. 2.2 in reproducing the code results in terms of integral quantities of importance, like the total number and total energy injected into the pair distribution.
The total number of Bethe-Heitler pairs produced per unit time, in the case of monoenergetic proton and photon populations, is given by dN BH,tot dt tot = N p N ph q BH (γ e ) dγ e (A.2) where N p and N ph represent the total number of protons and photons, respectively.In a similar manner the total amount of dimensionless energy that is injected into relativistic pairs is given by dE BH,tot dt tot = N p N ph q BH (γ e )γ e dγ e m e c 2 (A.3) Inspection of both panels in figure 20 shows that the empirical function of Eq. (2.2) reproduces satisfactorily the numerical results for the various test cases displayed.In particular, we recover the increasing trend in the number of injected pairs and energy with the proton Lorentz factor (yellow symbols and lines).On the contrary, when the proton Lorentz factor is fixed but the photon energy increases, both the number of pairs and the energy transferred to them drops as we move away from the threshold (see blue and red markers).Based on these results we can understand that the increasing trend found for the cases plotted in yellow is due to the presence of more energetic protons and not to the higher interaction energy, γ p ϵ, itself.This result might be explained from the fact that the energy transferred to the secondary population is a fraction of the proton energy [26], so as the latter increases, the energy transferred to the pair population will follow the same trend.
While the results depicted with yellow and blue symbols are obtained for the same protonphoton energies (see figure legend), the trend with interaction energy is different because of the number of interacting protons: the number of the injected protons remains the same for all cases shown in yellow, while for the cases shown in blue the compactness of the proton population ℓ p is kept fixed.Since N p,tot ∝ l p γ −1 p , the combination of a constant ℓ p with an increasing γ p means that the total number of protons decreases.On the contrary if N p,tot is constant and γ p increases, then ℓ p will also be increasing.Since the energy of the pair population is a fraction of the proton population energy, we expect that the energy transferred to pairs will be higher for cases with fixed N p,tot (yellow symbols).
Moreover it is interesting to compare cases of fixed proton Lorentz factor and variable target photon energy (red symbols) with cases where the proton Lorentz factor changes and the target- photon energy is constant (blue symbols).In both examples the compactnesses of both populations are fixed.We notice, that the total number of the produced pairs as a function of the interaction energy, γ p ϵ, is the same for both cases, see top panel in figure 20.However, the energy transferred to the pair population is significantly higher in the case where γ p increases, see blue symbols in the bottom panel of figure 20.This observation is expected if we consider that the total energy the pairs acquire (i) decreases away from the threshold, and (ii) is a fraction of the total energy of the proton population.In both the red and the blue cases, the interaction energy moves away from the threshold, thus we expect the total energy transferred to the pair population to decrease as γ p ϵ increases.This is, indeed, in the red cases where the photon-target energy is increased while the proton Lorentz factor remains constant.On the contrary, in the cases where γ p increases with ϵ being fixed, the total energy acquired by the pairs remains almost constant.
Lastly if we take a closer look at the magenta points we see that the total number of pairs is  approximately the same for interactions close to the threshold (γ p ϵ ≈ 2), even when γ p and ϵ are very different.However, the total energy of the pair population increases for higher proton Lorenz factors.This is in agreement with the fact that the pairs receive a fraction of the proton energy [26], and as the latter increases, the former will too.In other words, more energetic proton populations interacting with photon fields close to the threshold will produce approximately the same number of pairs, but with higher energy on average.

A.3 Comparison to other numerical implementations of the Bethe-Heitler injection spectrum
In this work we have constructed an analytical function that approximates the injection rate of the particle distribution created due to Bethe-Heitler pair production.The aforementioned function was benchmarked against the ATHEνA code [18,26].However, Kelner and Aharonian (KA) [12] provided with a formula that describes the Bethe-Heitler distribution of electrons or positrons produced when a single proton interacts with a photon population.This formula involves the numerical computation of three integrals, as shown in Eq. (2.1) [6].In this section, we use the leptohadronic numerical code LeHaMoC [21] that implements and solves the Kelner and Aharonian formula for Bethe-Heitler pair production, in order to compare our analytical function against the aforesaid formula.
Figure 21 shows the comparison of the injection spectrum obtained with the analytical function we have introduced in this work (Eq.2.2), the results of the benchmarking code ATHEνA, and the results of the KA formula as implemented in LeHaMoC.We perform this comparison for two different parameter sets.First, we tested the interaction of a monoenergetic proton population for three different proton Lorentz factor values, γ p = 10 4 , 10 6 , 10 7 , with a grey body target-photon distribution of temperature T = 10 6 K (left panel).Second, in the right panel, we show the comparison of the three methods for a power-law proton distribution, n p = n p,0 γ −2 p with γ p ∈ [10 4 , 10 8 ] interacting with a power-law target photon field, n ph = n ph,0 ϵ −2 with ϵ ∈ [10 −8 , 10 −4 ].For all the curves we use a proton compactness of l p = 10 −4 and a photon compactness of l ph = 8 • 10 −6 for the grey body photons and l ph = 10 −1 for the power-law photon distribution.Our analytical approximation matches satisfactorily the ATHEνA results as expected (see previous sections).More importantly, one notices that the analytical function is in a good agreement with the KA formula in most cases (blue and yellow lines of the left panel and magenta line of the right panel).However, for the red lines of the left panel of figure 21, where γ p = 10 7 , we see that the KA formula under predicts the production rate of pairs compared to ATHEνA, and as a result does not match the results of our analytical function.The latter is due to the fact that the KA formula is not a valid description of the Bethe-Heitler created pair population when the interaction energy, γ p ϵ, reaches values comparable to the mass ratio of protons and electrons, m p /m e , which is the case for the red lines 5 of the left panel in figure 21.This mismatch does not play in important part in the case of extended proton and/or photon distributions (i.e.right panel of figure 21) because in such cases the characteristics of the pair distribution are determined of the near-threshold interactions for which the three lines are in good agreement (blue curve of the left panel).

B Synchrotron cooling for different B values
In section 2, we presented steady-state cooled distributions of Bethe-Heitler produced pair populations in a source having a magnetic field strength B = 40 G (see figure 4).In figure 22 we present steady-state pair distributions with the synchrotron spectra each population emits, for different cooling regimes.
Figure 22 shows that in sources with high magnetic field values, B > 4 G, the maximum of the synchrotron spectra does not appear at Lorentz factor value γ e ≈ γ br .On the contrary, both the pair energy distribution and the electromagnetic spectrum peak at either γ e ≈ γ p,min or γ e ≈ γ p,max , depending on the slope of the target photon field (see figure 3).However, in sources with weaker magnetic fields, B < 4 G, we observe that γ e,min ≤ γ br ≤ γ e,max , so the particles with γ e < γ br do not cool significantly in a time frame of R ′ b /c, where R ′ b = 10 15 cm in these examples.As a result the steady-state energy distribution is a combination of the injected energy distribution, as shown in figure 2, with an extra peak originating by the concentration of cooled high-energy pairs around γ e ≈ γ br .The latter feature is more visible for photon power-law slopes s ph ≤ 1 and it is, also reflected, in the synchrotron spectra of the population.

C Additional blazar SEDs
In section 3.3 we presented results for a specific source radius.Blazar SEDs computed for smaller and larger emitting regions are presented in figures 23 and 24.For larger emitting volumes (∝ R ′,3 b ) higher photon number densities are required for making the source optically thick to γγ pair production.Thus, even by keeping the proton and electron injected luminosities values fixed, we can achieve higher photon luminosity values by increasing the size of the blob.In figures 23 and 24 we observe that there are no cases that close to the γγ optically thick source limit.When it comes to the variation of ℓ e and ℓ p , the behavior of the large sources is similar to the one found for R ′ b = 10 15 cm.

D SED model parameters
Input parameters to the ATHEνA code describing the injection populations of electrons, protons and, in some cases, external photons leading to the theoretical models of section 3.

Figure 1 :
Figure 1: Injection spectra, produced by ATHEνA [18, 26], of Bethe-Heitler pairs produced by interactions of (i) a single proton with γ p = 10 5.5 with a single photon of different energies (see color bar in panels a and c) and (ii) a single proton of a changing Lorentz factor (see color bar in panel b) with a single photon of energy adjusted, in each case, in order to keep the interaction energy constant.Panel (d) shows the maximum value of the energy injection spectra from panels (b) and (c) as a function of the proton Lorentz factor (normalized to 10 5 ) and the interaction energy, respectively.

Figure 2 :
Figure 2: Injection spectra of pairs produced by interactions of power-law protons with slope s p = 2 and γ p ∈ [10 4 , 10 8 ] interacting with power-law photons with ϵ ∈ [10 −8 , 10 −4 ] and different photon indices (see color bar).Solid lines represent numerical results from the ATHEνA code [18, 26] and dashed lines represent results obtained with the analytical approximation given by Eq. (2.2).To facilitate the comparison of the two solutions, a ratio plot is added below the main panels, with the lower panel zooming into values in the range of 0.8-1.2 .

Figure 3 :
Figure 3: Energy injection spectra of pairs produced in Bethe-Heitler interactions between protons with a power-law (PL) distribution for different slopes s p , and PL photon distributions with different photon indices, s ph .All other parameters are the same as in figure 2. Black solid lines represent the total injection spectrum obtained after convolving the empirical function of Eq. (2.2) with the proton and photon distributions.Colored solid lines represent the contribution of single-energy protons (see color bar) to the total spectrum.When s ph − s p > 1 (< 1) the peak position and value of the energy injection spectrum is determined by the highest (lowest) energy protons of the distribution.

Figure 5 :
Figure 5: Analytical constraints for γ-ray production by Bethe-Heitler pairs following from Eq. (3.4) (shaded regions) and Eq.(3.5) (colored lines).Protons interact at the threshold with electron synchrotron photons, and produce pairs that emit synchrotron radiation within the observed γ-ray limits (see Table3) for γ ′ e and γ ′ p values selected from the intersection of colored lines with the shaded region for each blazar subclass.The peak energy of the SSC component is lower than the peak energy of the Bethe-Heitler synchrotron component for parameter values drawn from the magenta shaded region.There are combinations of γ ′ p and γ ′ e leading to Bethe-Heitler synchrotron spectra peaking in γ rays.

. 11 )
leaving more freedom in the determination of γ ′ e and γ ′ p .If we, now, combine Eqs.(3.1) and (3.11), we can express the temperature of the external thermal photon field as:

Figure 6 :
Figure 6: Graphical representation of the γγ optical depth as a function of the high-energy photon energy for different blazar classes.Colored areas indicate the typical range of observed high-energy peak energies.For each class we integrate Eq. (3.27) on an indicative blazar spectrum.We use the photon spectra of figure 8 for B ′ = 50 G.The latter are also shown in the figure with dashed lines, in the blob frame and are normalized to their peak values.We also take into account the Doppler factor and the redshift of each source, as δ = 55, 45, 27 and z = 0.5, 0.3, 0.1 for the LSP, the ISP and the HSP cases respectively.In all cases, the emitting region is opaque to γγ pair production for photons of energy above 10 TeV.

Figure 7 :
Figure 7: Decomposition of an indicative spectral energy distribution (SED) from an ISP-like blazar into various spectral components indicated by different types of lines.The all-flavor neutrino spectrum is also plotted (black solid line).Colored markers on the upper horizontal axis indicate analytical estimates of characteristic energies presented in section 3.1.For the parameter values used, see Table 4 under ISPs with R ′ b = 10 14 cm and B ′ = 50 G.

Figure 8 :
Figure 8: Theoretical spectral energy distributions (SEDs) for LSP-, ISP-, and HSP-like blazars (thick solid lines) computed for an emitting region with radius R ′ b = 10 15 cm and different values of the (comoving) magnetic field strength (all other parameter values are listed in Table4).In all panels, dashed lines represent the synchrotron spectrum of Bethe-Heitler pairs computed using the analytical function for the pair injection, see Eq. (2.2).Dash-dotted and dotted lines represent spectra obtained by varying the proton and electron compactness, respectively, by two orders of magnitude with respect to their nominal values (solid lines).The energies of the low-and high-energy humps for the baseline model (solid lines) are indicated with symbols.Long-dashed lines in the LSP graphs represent the external photon field.Gray shaded regions indicate the observational range of peak frequencies for the low-and high-energy SED components.Black horizontal lines at ∼ 1 − 10 PeV represent the neutrino flux prediction according to Eq. (3.26).

Figure 9 :
Figure 9: Bolometric luminosities of the primary electrons and Bethe-Heitler injected pairs and their emitted synchrotron radiation.Panel (a) shows the ratio of the synchrotron luminosity to the injected luminosity derived numerically for primary electrons (filled symbols) and secondary pairs (open symbols).Panel (b) displays the ratio of the primary electron synchrotron luminosity computed numerically to the analytical estimate of using Eq.(3.16).Panel (c) shows the same ratio as in panel (b) but for Bethe-Heitler injected pairs, with the analytical estimate given by Eq. (3.20).All the ratios are plotted against the synchrotron cooling timescale of particles injected with Lorentz factor γ e to the timescale R ′ b /c.

Figure 10 :
Figure 10: Peak synchrotron luminosity (left panel) and Compton dominance (right panel) versus the peak synchrotron frequency.Our model results are plotted with colored markers (see inset legends) on top of observationally inferred values (gray crosses) from[51].The size of colored markers indicates the size of the blob, with the smallest markers corresponding to 10 14 cm and the largest to 10 16 cm.Open markers indicate models that (i) fall outside the locus of observed points, or (ii) have γ-ray peak frequencies outside the typical range of values for each subclass (see table3), or (iii) require non-physical combinations of temperature and luminosity of the external photon field (see section 3.3.3).

Figure 11 :
Figure 11: Chart of number of total models (black) and valid models (red), categorized in respect of their electromagnetic spectrum type (LSP-, ISP-, HSP-like).

Figure 12 :
Figure 12: Jet power (Eq.3.29) as a function of the ratio of the relativistic proton and magnetic field energy densities (left panel), and the peak γ-ray luminosity of the source (right panel).Smaller size and more transparent markers indicate smaller emitting regions.The horizontal line marks the Eddington luminosity of a 10 8 M ⊙ black hole.

Figure 13 :
Figure 13: Graphical representation of the three exponential terms appearing in Eq. (2.2) without accounting for the normalization A. All curves in the left panel are normalized to a peak value of one.The magenta line in both panels represents the product of the three terms.

Figure 14 :
Figure 14: Left panel: Power-law index p appearing in Eq. (2.2) as a function of the interaction energy γ p ϵ.The dashed blue line represents the MCMC fit result with the gray shaded area being the 1σ spread of the values.The solid magenta line is another solution obtained for parameter values slightly different than the median values of the posterior distributions (see corner plot on the right panel).Right panel: Corner plot showing the posterior distributions of the parameters in Eq. (A.1).Dashed lines in the histograms indicate the median value and the 68% confidence interval.Solid black lines represent the values that reproduce the magenta line of the left panel.

Figure 15 :
Figure 15: Energy injection spectra for the cases shown in figure 1. Solid and dashed lines show results when p(γ p ϵ) follows the the magenta and blue curves, respectively, shown in figure 14.Greyshaded areas demonstrate spectra obtained for values of p(γ p ϵ) drawn from the 68% of figure 14.

Figure 16 :
Figure 16: Normalization A of the injection function (see Eq. 2.2) plotted against the proton Lorentz factor γ p , or the target-photon energy ϵ, or the interaction energy γ p ϵ. Red symbols show results when the proton Lorentz factor changes but the photon energy is ϵ ≈ 2/γ p (i.e.interactions happen close to the threshold).Magenta symbols show results when photons of different energies interact with protons of a fixed Lorentz factor (i.e. the interaction energy changes).Blue symbols show results when protons of different Lorentz factors interact with photons of ϵ = 10 −5 .In all cases, the proton and photon distributions are normalized to their total number.The grey shaded areas represent the normalization function range in the 1σ spread of the fitted values.

Figure 17 :
Figure 17: Corner plot showing the posterior distributions of the parameters of the normalization function A displayed in Eq. (2.6).Dashed lines in the histograms indicate the median value and the 68% confidence interval.Solid black lines show the values reported in table 2.

Figure 18 :
Figure 18: Injection spectra of Bethe-Heitler pairs produced by interactions of a single proton with γ p ≈ 3 • 10 5 with a single photon of different energies (see color bar) -same cases as in panel c of figure 1. Solid lines represent numerical results from the ATHEνA code [18, 26] and dashed lines represent the empirical Bethe-Heitler spectrum of Eq. (2.2).Bottom panels show the ratio of the code to analytical results.

Figure 19 :
Figure 19: Same as in figure 18 but for different proton Lorentz factors (see color bar).The target photon energy in each case is adjusted in order to keep the interaction energy constant, i.e. γ p ϵ ≃ 5.

Figure 20 :
Figure 20: Integrals of total number (upper left panel) and energy (upper right panel) of injected Bethe-Heitler pairs per unit time (with residuals comparing against the ATHEνA results) in cases where (i) the proton Lorentz factor changes while the photon energy is ϵ ≈ 2/γ p ; larger marker size indicates higher γ p values (magenta), (ii) photons of different energies interact with protons of the same energy (red), (iii) protons of different Lorentz factors interact with 10 −5 energy photons when the former compactness changes so the total number of the injected protons to be constant (blue) or l p ∝ γ p (yellow).

Figure 22 :
Figure 22: Steady-state energy distributions of Bethe-Heitler pairs and their emitted synchrotron spectra for various magnetic field strengths and different power-law slopes of the target photon field (see color bar).Vertical dotted lines show the synchrotron emission of the pairs with the lowest and highest Lorentz factor values, while vertical magenta lines represent the synchrotron cooling break Lorentz factor, defined in Eq. 2.9.
To ensure continuity of q BH at γ e = 15γ p we multiply Eq. (2.2) with

Table 4 :
6.5 ,7.0] Parameter values for the injected populations and the source characteristics for the all cases