Bethe-Heitler Cascades and Hard Gamma-Ray Spectra in Flaring TeV Blazars: 1ES 0414+009 and 1ES 1959+650

1. Introduction

Active galactic nuclei (AGNs) are among the most numerous extragalactic objects in astronomy, with blazars representing an extreme subclass of radio-loud AGNs. Blazars exhibit non-thermal continuum emission originating from jetted AGNs, with the jet axis closely aligned to the observer’s line of sight [1]. Moreover, they are characterized by distinctive characteristics such as rapid variability, high luminosity, and strong polarization [2,3,4,5,6]. Blazars are classified into two subclasses: flat-spectrum radio quasars (FSRQs) and BL Lac objects. FSRQs are distinguished by their strong, broad emission lines, whereas BL Lac objects exhibit weak or nearly absent emission lines [7]. Their multiwavelength spectral energy distributions (SEDs) typically feature a characteristic double-hump structure [3,4,8,9,10,11]. The low-energy hump, spanning from the radio to the X-ray bands, is well explained by synchrotron radiation. However, the origin of the high-energy hump, which falls within the MeV-TeV energy range, remains an open question. Two theoretical models describe the high-energy photon emission in blazars: the leptonic model and the hadronic model. In the leptonic model scenarios, the high-energy component may result from IC scattering of relativistic electrons, either on synchrotron photons (Synchrotron Self-Compton, SSC [12,13]) or on external photon fields (External Compton, EC [4,14,15]), such as the broad-line region (BLR), accretion disk, or the external torus. In contrast, the hadronic model suggests that high-energy gamma-rays are generated either through proton synchrotron radiation in sufficiently strong magnetic fields [16,17,18] or via meson and lepton production in cascades triggered by proton-proton or proton-photon interactions [19,20,21,22].

The position of the first peak in SEDs, ${\nu }_p^S$ (synchrotron peak frequency), classifies the sources as low-frequency peaked BL Lacs (LBLs; e.g., ${\nu }_p^S$ < ${10}^{14}$ Hz), intermediate-frequency peaked BL Lacs (IBLs; e.g., ${10}^{14}$ Hz <${\nu }_p^S$ < ${10}^{15}$ Hz), and high-frequency peaked BL Lacs (HBLs; e.g., ${\nu }_p^S$ > ${10}^{15}$ Hz) [8]. Several studies [12,23,24] have shown that the SEDs of BL Lac objects, particularly HBLs, are well described by a pure SSC model. The formation of relativistic jets in AGNs remains an open question, with various models proposed to explain their origin. Two of the most well-established theories are the Blandford-Znajek mechanism [25], in which the jet derives energy from the black hole’s rotation, and the Blandford-Payne mechanism [26], where the jet primarily extracts rotational energy from the accretion disk. In both scenarios, magnetic fields play a crucial role in channeling energy from the black hole or disk into the jet [27]. Initially, the jet’s energy is primarily carried by Poynting flux, which progressively converts into the kinetic energy of the plasma as the flow accelerates [28]. Theoretical studies on energy dissipation suggest that, in this scenario, the emitting electrons and the electromagnetic field may share an equal distribution of energy [29]. Understanding the strength of the magnetic field within the jet is essential for unraveling its formation and energy distribution. By analyzing the frequency-dependent position of the optically thick jet core [30], Zamaninasab et al. [31] found that the jet’s magnetic flux on parsec scales correlates with the power of the corresponding accretion current, aligning with predictions from magnetohydrodynamics. Another reliable approach for estimating the magnetic field in the innermost emission region of the jet is through modeling the SEDs of blazars [3,32,33,34].

Comparing multi-wavelength emission from blazars with numerical models that simulate radiative emission and transfer under specific assumptions about particle content and emission region characteristics is a key method for investigating the source’s microphysics and physical parameters. Data from the very high-energy (VHE) segment of the SED is essential for constraining model parameters of high-frequency peaked BL Lac objects, which emit a significant portion of their radiation in this energy range. Additionally, VHE observations are valuable for softer sources, as they probe the extreme limits of particle acceleration and are particularly sensitive to photon-photon absorption within internal or external radiation fields [35]. An essential element for understanding the high-energy gamma-ray emission is the Bethe-Heitler pair production (also called the proton-photon ($p\gamma$) pair production process), which occurs when relativistic protons interact with soft photons (ambient synchrotron or external photons), leading to the production of electron-positron pairs [36,37]. The generated electron-positron pairs then undergo synchrotron and IC processes, contributing to the broadband gamma-ray emission observed from blazars. It is worth mentioning that this process is typically subdominant compared to proton synchrotron or pion decay (which predict correlated neutrino and gamma-ray emissions) mechanisms in hadronic models but can still play a role in shaping the observed SED of blazars.

In this work, we investigate the SED of two distant blazars detected at TeV energies: 1ES 0414+009 and 1ES 1959+650. These high-energy sources provide valuable insights into particle acceleration mechanisms and emission processes. By modeling their multiwavelength emission, we aim to constrain the physical properties of their jets and explore the contribution of leptonic and hadronic processes to their observed radiation. In Section 2, we introduce our samples. In Section 3 and Section 4, we present the fitting tools and the fitting process. Section 5 describes the results and provides a discussion. Finally, we present the conclusion of this work in Section 6. Throughout the paper, we assume the Hubble constant $H_0$ = 75 km ${\mathrm{s}}^{-1}$ ${\mathrm{Mpc}}^{-1}$ , the matter energy density ${\Omega }_M$ = 0.27, the radiation energy density ${\Omega }_r$ = 0, and the dimensionless cosmological constant ${\Omega }_{\Lambda }$ = 0.73.

2. Distant TeV Blazars: 1ES 0414+009 and 1ES 1959+650

In this section, we are going to present a brief description of the two distant TeV blazars studied in this work and the set of data used in our analysis.

2.1. 1ES 0414+009

The BL Lac object 1ES 0414+009 was first detected by the HEAO 1 satellite [38] in the 0.2 keV–10 MeV energy range and later identified in X-ray images from the Einstein Observatory [39]. Situated at a redshift of z = 0.287 [40], it is powered by a supermassive black hole (SMBH) with a mass of approximately 2 × ${10}^9$ ${\mathrm{M}}_{\odot }$[41]. According to the scheme proposed by [42], 1ES 0414+009 belongs to the class of HBLs, which exhibit a synchrotron-emission peak at UV/soft X-ray frequencies. In such sources, X-ray emission is primarily dominated by synchrotron radiation. In the VHE gamma-ray domain, data from the HEGRA experiment were used to establish an upper limit on the flux of 1ES 0414+009, corresponding to 13.5 × ${10}^{-12}$ ${\mathrm{cm}}^{-2}$ ${\mathrm{s}}^{-1}$ above 910 GeV [43]. Costamante & Ghisellini [44] identified 1ES 0414+009 as a strong candidate for VHE emission based on its high X-ray and radio flux. Its detection became even more probable following blazar gamma-ray spectrum analyses, which suggested a low intensity of the diffuse extragalactic background light (EBL) [45].

The H.E.S.S. array of Cherenkov telescopes has detected significant VHE gamma-ray emission from 1ES 0414+009. With an average flux of ∼0.6% of the Crab Nebula flux above E > 200 GeV, this blazar is among the faintest extragalactic sources observed in the TeV range. Additionally, 1ES 0414+009 was detected by the Fermi-LAT instrument during its first 20 months of operation (2008–2010), exhibiting very faint emission in the high-energy (HE) domain as well [46]. The HE and VHE spectra of 1ES 0414+009, corrected for absorption using an EBL model close to the lower limits, exhibit a best-fit power law with an index harder than 2, classifying it as a hard-TeV BL Lac object. The overall SED is averaged over five years, though there is no strict simultaneity between Swift and H.E.S.S. observations. With this limitation, the SED properties, particularly an IC peak energy above 1–2 TeV, are challenging to explain within the framework of a pure one-zone SSC model, unless unusual parameter values are assumed [47].

The study of 1ES 0414+009 has implications beyond its characterization, extending to the broader context of AGN and their role in multi-messenger astronomy [48]. This blazar, along with others, has been proposed as a potential source of ultra-high-energy cosmic rays [49]. These investigations collectively enhance our understanding of the extreme environments surrounding supermassive black holes and their relativistic jets, contributing to advancements in high-energy astrophysics [50,51].

2.2. 1ES 1959+650

The object 1ES 1959+650 [52], with a redshift of z = 0.047, was classified as a BL Lac in 1993 using a specialized radio/optical/X-ray technique [53]. Its first detection at TeV energies was reported by the Utah Seven Telescope Array collaboration during the 1998 observational season [54], revealing an excess with a statistical significance of 3.9$\sigma$ above 600 GeV after 57 hours of observation. In May 2002, 1ES 1959+650 experienced a powerful TeV outburst, observed by the VERITAS [55], HEGRA [56] collaborations. Significant flux variations were recorded, reaching levels up to three times the Crab Nebula flux [57].

As a HBL source, 1ES 1959+650 is generally faint in the Fermi-LAT energy range (20 MeV–300 GeV) compared to nearby low-energy-peaked BL Lacs (LBLs) and exhibits weaker variability in this band than in the X-ray and VHE ranges [58]. Notably, the 0.3–10 GeV flux did not display significant long-term flares, with the photon flux derived from two week binned data rarely exceeding 4 × ${10}^{-8}$ photons ${\mathrm{cm}}^{-2}$ ${\mathrm{s}}^{-1}$ [59,60,61]. However, the source remained mostly above this threshold between August 2015 and August 2016, during which two strong, long-term HE flares were observed [62,63].

During intense X-ray flares in 2016–2017, the source exhibited very hard X-ray spectra, with the 0.3–300 GeV photon index also remaining hard during the same periods [63]. Notably, achieving a hard gamma-ray spectrum is more naturally explained within hadronic scenarios, whereas reproducing such a spectrum with leptonic models is significantly more challenging [64]. For instance, the proton blazar model, introduced by Mannheim [20], predicts X-ray spectra with a photon index of 1.5–1.7 and an uncorrelated X-ray–TeV variability, a characteristic observed in our target, which is more easily accounted for by hadronic models [58].

3. Multiwavelength SED Fitting of Jetted AGNs

Modeling the spectral energy distributions (SEDs) of jetted active galactic nuclei (AGNs) across the electromagnetic spectrum, from radio to gamma rays, is essential for understanding the physical mechanisms governing relativistic jets. Through multiwavelength SED fitting, we can investigate particle acceleration, energy dissipation, and radiative processes responsible for the observed emission.

3.1. SED Leptonic Modelling

We employed the JetSeT software package [65,66,67] for leptonic SED modeling. JetSeT is an open-source C/Python framework designed to simulate radiative and particle acceleration processes in relativistic jets and galactic sources, both beamed and unbeamed. It supports various leptonic emission mechanisms, including synchrotron radiation, synchrotron self-Compton (SSC), and external Compton (EC) scattering of photons from the accretion disk, broad-line region (BLR), dusty torus (DT), and the cosmic microwave background (CMB). The framework also includes $\gamma \gamma$ absorption based on established extragalactic background light (EBL) models [68,69,70].

To ensure consistency across sources, we adopted a one-zone synchrotron+SSC scenario, where the emission originates from a spherical region of radius $R=c{t}_{\mathrm{var}}\delta /(1+z)$ [71], with $t_{\mathrm{var}}$ set to one day [72,73]. The region moves relativistically with Doppler factor $\delta ={[\Gamma (1-\beta cos\theta )]}^{-1}$, where $\Gamma$ is the bulk Lorentz factor and $\theta$ the viewing angle. Electrons are assumed to follow a broken power-law distribution [12,23,28]:

$$N\left(\gamma \right)=N_0\left\{\begin{array}{cc}{\gamma }^{-p_1},\hfill & {\gamma }_{\min }\le \gamma \le {\gamma }_{\mathrm{b}},\hfill \\ {\gamma }_b^{p_2-p_1}{\gamma }^{-p_2},\hfill & {\gamma }_{\mathrm{b}}<\gamma \le {\gamma }_{\max },\hfill \end{array}\right.$$

(1)

where ${\gamma }_{\min }$, ${\gamma }_{\mathrm{b}}$, and ${\gamma }_{\max }$ are the minimum, break, and maximum electron Lorentz factors, and $N_0$ is the normalization constant.

JetSeT employs a two-stage fitting procedure. First, a phenomenological characterization of the SED is performed using the SEDShape module, which applies power-law and log-parabolic fits to binned data spanning the radio to TeV range. This step provides initial constraints on the synchrotron and SSC components, helping to define parameter boundaries. Next, the ObsConstrain module is used to derive input parameters for the physical modeling, where we adopt the broken power-law electron distribution defined in Equation (1). A Bayesian approach is adopted for the final model fitting, incorporating prior constraints to ensure that all parameters remain within physically meaningful bounds.

Model parameters are implemented as Astropy quantities, allowing seamless integration with other Python-based astrophysical tools. JetSeT supports both frequentist and Bayesian fitting approaches. For frequentist optimization, plugins are available for iminuit [74] and SciPy’s bounded least squares method [75]. Bayesian inference is carried out using an MCMC sampler through integration with the emcee package [76]. In all cases, the redshift is fixed, and we assume a cold proton to relativistic electron ratio of 0.1 [77].

3.2. SED Lepto-Hadronic Modelling

To explore scenarios involving hadronic contributions, we used the ${\mathrm{AM}}^3$ (Astrophysical Multi-Messenger Modeling) framework [78], an open-source software package designed to simulate time-dependent lepto-hadronic interactions in astrophysical sources. ${\mathrm{AM}}^3$ computes the coupled evolution of photon, electron, positron, proton, neutron, and neutrino populations, along with intermediate products, within an isotropic magnetic field. It includes non-linear processes such as electromagnetic cascades and secondary particle feedback, providing a self-consistent description of particle interactions.

The emission region is modeled as a spherical blob of radius $R^{\prime }$ moving with bulk Lorentz factor $\Gamma$ along the jet. Primary electrons and protons are injected isotropically into this region. Electrons follow a broken power-law energy distribution, while protons follow a single power-law extending up to ${\gamma }_{p,\max }^{\prime }$. These high-energy protons interact with ambient photon fields following the framework of Hümmer et al. [79], producing charged and neutral pions. The decay of pions gives rise to secondary gamma rays, neutrinos, electrons, and positrons, which in turn participate in electromagnetic cascades.

${\mathrm{AM}}^3$ incorporates several key processes: synchrotron emission and self-absorption, IC scattering by both electrons and protons, Bethe–Heitler pair production ($p+\gamma \to p+e^{-}+e^{+}$), photon-photon pair production and annihilation, and the evolution of secondary particles. The magnetic field $B^{\prime }$ is assumed to be randomly oriented within the blob, and its turbulence plays a significant role in shaping the resulting SED and multi-messenger signatures.

This approach allows us to assess the hadronic contribution to the SED and explore scenarios where neutrino and gamma-ray emission arise from the same physical region, providing insight into possible associations between high-energy astrophysical neutrinos and blazar flares.

4. Results and Discussions

Based on the methods presented in Section 3.1 and Section 3.2, it was possible to perform the non-thermal multiwavelength modeling of the blazars 1ES 0414+009 and 1ES 1959+650. For this, data from space and ground-based observatories, ranging from radio to VHE gamma-rays, were extracted from the Space Science Data Center (SSDC)1 and Firmamento2. Firstly, the multiwavelength observations were fitted using the MCMC method in JetSeT, and after the best-fitting parameters were used as input in the lepto-hadronic modeling via ${\mathrm{AM}}^3$. Table 1 and Table 2 summarize the observatories employed in the modeling of the sources 1ES 0414+009 and 1ES 1959+650, respectively, together with the corresponding energy ranges covered by the data collected at each facility.

4.1. SSC SED Fitting

In the fitting method, each free parameter has a specific physical boundary. We can fix specific parameters and set the fitting range of the remaining parameters to speed up the convergence of the fitting process. Given that all our samples are HBLs and we adopt the one-zone lepton model Syn+SSC, we fix the redshift and the distance of the radiation region from the central black hole, $R_H$ - ${10}^{17}$ cm (default value for JetSeT). To prevent biased results, ensure that output parameters remain within physically meaningful ranges (based on the literature [81] and [107]), and improve convergence, we define fitting ranges for the radius of the emission region R, magnetic field strength B, ${\gamma }_{\min }$, ${\gamma }_{\mathrm{b}}$, ${\gamma }_{\max }$, and the spectral indexes (p and $p_1$). The best-fitting values are described in Table 3.

Figure 1 shows the modeling via the SSC physical process for the 1ES 0414+009 source. The first peak in the figure represents the synchrotron emission, which occurs at lower frequencies, peaking at $\sim {10}^{17}$ Hz, thus indicating that 1ES 0414+009 is a HSP BL Lac. The second peak, at $\sim {10}^{25}$ Hz, represents the emission due to IC scattering. According to the residual plot at the bottom of the figure, the model fits the data well by plotting the difference between actual and predicted values, with only one discrepancy in the TeV region of the energy spectrum.

In addition to this, the blazar jet contains both the energy density ($U_e$) of relativistic electrons and the energy density of the magnetic field ($U_B$) at the blob rest frame. For 1ES 0414+009, under the simple one-zone lepton Sync+SSC model, $U_e\approx 3.22\times {10}^{-6}$ erg/${\mathrm{cm}}^3$ and $U_B\approx 1.60\times {10}^{-4}$ erg/${\mathrm{cm}}^3$. Moreover, the radiation energy density of synchrotron photons is $U_{\mathrm{Synch}}\approx 2.169\times {10}^{-8}$ erg/${\mathrm{cm}}^3$. In JetSeT, jet power can be readily estimated by fitting the spectral energy distributions. The jet kinetic power is carried by electrons, magnetic fields, and cold protons [28,71,108], and is expressed as $L_i=\pi R^2{\Gamma }^{2}c{U}_i$, where $U_i$ represents the energy density of each component ($i=e,p,B$). Consequently, the total jet kinetic power is given by $L_{\mathrm{kin}}=L_e+L_B+L_p$. The derived energy densities and luminosities for the source 1ES 0414+009 are summarized in Table 4. The total radiative power is calculated as $L_{\mathrm{rad}}\simeq L_{\mathrm{syn}}+L_{\mathrm{SSC}}$, where $L_{\mathrm{syn}}$ and $L_{\mathrm{SSC}}$ are the powers emitted through synchrotron and synchrotron self-Compton processes, respectively [13]. For 1ES 0414+009, we find $L_{\mathrm{syn}}\approx 1.62\times {10}^{39}$ erg ${\mathrm{s}}^{-1}$ and $L_{\mathrm{SSC}}\approx 5.18\times {10}^{35}$ erg ${\mathrm{s}}^{-1}$, resulting in a total radiated power of $L_{\mathrm{rad}}\simeq 1.62\times {10}^{39}$ erg ${\mathrm{s}}^{-1}$. Based on the energy densities, the jet kinetic power is estimated as $L_{\mathrm{kin}}=(2.44\times {10}^{41}+1.21\times {10}^{43}+2.82\times {10}^{43})\approx 4.06\times {10}^{43}$ erg ${\mathrm{s}}^{-1}$.

The comparison between the parameters of the SSC model for 1ES 0414 + 009 obtained in this work and those reported by Aliu et al. (2012) [109] reveals important similarities despite differences in the modeling approaches and data selections. Both studies find a comparable size for the emitting region, with $R\sim {10}^{16-17}$ cm, consistent with a compact emission zone typical of high-frequency peaked BL Lac objects. In addition, the viewing angle (${\theta }_{\mathrm{obs}}\sim 1.4^{\circ }$) and the general SSC framework adopted for the SED modeling are similar, indicating a shared assumption of a relativistic jet closely aligned with the observer’s line of sight. Both works also conclude that purely leptonic models face challenges in fully reproducing the broadband SED, particularly in explaining the hard TeV spectra, thus motivating the consideration of more complex scenarios, such as lepto-hadronic contributions. Although there are differences in specific parameter values, such as the minimum Lorentz factor, magnetic field strength, and bulk Lorentz factor, these arise naturally from the focus of each study: Aliu et al. emphasized modeling the TeV emission with extreme parameters, while our work provides a global fit to a broader, multiwavelength dataset. Overall, both studies reinforce the notion that 1ES 0414+009 exhibits characteristics requiring detailed modeling beyond the simplest SSC assumptions.

Similarly to the previous source, Figure 2 shows the modeling of 1ES 1959+650 using the SSC physical process, considering also the presence of the host galaxy component (emission peak at approximately ${10}^{15}$ Hz). As indicated by the figure, the synchrotron emission peaks at $\sim {10}^{18}$ Hz, classifying the source as an extreme high-frequency-peaked BL Lac object (EHBL) with ${\nu }_{\mathrm{sync}}^{\mathrm{peak}}>{10}^{16}$ Hz. The second peak, located at higher energies, is attributed to IC scattering. For 1ES 1959+650, the best-fitting parameters are ${\gamma }_{\min }=2.78\times {10}^2$, ${\gamma }_{\mathrm{break}}=1.019\times {10}^3$, ${\gamma }_{\max }=1.34\times {10}^6$, $B=5.08\times {10}^{-2}$ G, and $R=7.20\times {10}^{15}$ cm, with a relativistic particle number density of $N=38.79$ ${\mathrm{cm}}^{-3}$. The low- and high-energy spectral indices are $p=1.670$ and $p_1=2.576$, respectively, as summarized in Table 3.

At the blob, the electron energy density is $U_e\approx 5.04\times {10}^{-4}$ erg ${\mathrm{cm}}^{-3}$, the magnetic energy density is $U_B\approx 1.03\times {10}^{-4}$ erg ${\mathrm{cm}}^{-3}$, and the synchrotron photon energy density is $U_{\mathrm{sync}}\approx 1.71\times {10}^{-7}$ erg ${\mathrm{cm}}^{-3}$. The corresponding synchrotron and SSC radiative powers are $L_{\mathrm{syn}}\approx 8.30\times {10}^{37}$ erg ${\mathrm{s}}^{-1}$ and $L_{\mathrm{SSC}}\approx 3.34\times {10}^{35}$ erg ${\mathrm{s}}^{-1}$, respectively, resulting in a total radiated power of $L_{\mathrm{rad}}\simeq 8.33\times {10}^{37}$ erg ${\mathrm{s}}^{-1}$. These quantities are listed in Table 4. The jet kinetic power is estimated as $L_{\mathrm{kin}}=(2.45\times {10}^{41}+4.98\times {10}^{40}+2.83\times {10}^{43})\approx 2.86\times {10}^{43}$ erg ${\mathrm{s}}^{-1}$. As found for most HBLs, the kinetic powers carried by electrons and protons ($L_e$ and $L_p$) dominate over the magnetic contribution $L_B$ [28].

A comparison between Table 3 and Table 4 from this work and Table 1 from Aliu et al. (2013) [59] highlights similarities in the physical parameters derived for 1ES 1959+650. Both studies model the broadband emission using a one-zone leptonic scenario dominated by synchrotron and SSC processes, adopting similar assumptions for the jet composition and emission region structure. The magnetic field strength ($B\sim 0.02$–$0.06$ G) and the size of the emission region ($R\sim {10}^{16}$–${10}^{17}$ cm) are consistent across the studies, suggesting comparable estimates for the energy balance within the jet. Additionally, both works report high minimum Lorentz factors for the relativistic electrons (${\gamma }_{\min }\sim {10}^2$–${10}^4$) and relatively hard injected electron spectral indices ($p\sim 1.7$–$2.0$), indicative of efficient particle acceleration. The inferred jet kinetic powers are also similar, with both analyses concluding that the power carried by particles (electrons and protons) largely exceeds the magnetic contribution, a common characteristic of high-frequency-peaked BL Lac objects. These similarities reinforce the model of the one-zone SSC interpretation for the quiescent state of 1ES 1959+650, as presented in both studies.

4.2. Lepto-Hadronic SED Modelling

The lepto-hadronic SED modelling was performed using the open-source ${\mathrm{AM}}^3$ software. The non-thermal lepton distribution in the source came from the JetSeT fitting process (section above), where the electron luminosity in the jet is $L_e\sim 2.44\times {10}^{41}$ erg/s for 1ES 0414+009 and $L_e\sim 2.45\times {10}^{41}$ erg/s for 1ES 1959+650. For the hadronic interaction process, the proton luminosity $L_p$ was estimated to be less than or equal to the Eddington luminosity ($L_p\le L_{\mathrm{Edd}}$) [110], where $L_{\mathrm{Edd}}=1.3\times {10}^{38}(M_{\mathrm{BH}}/M_{\odot })$ erg/s. Assuming that $M_{\mathrm{BH}}\sim 2\times {10}^9{\mathrm{M}}_{\odot }$ for 1ES 0414+009 [82], $L_{\mathrm{Edd}}\sim 2.52\times {10}^{47}$ erg/s. For 1ES 1959+650 we also assumed $L_p\le L_{\mathrm{Edd}}$ [110], where $M_{\mathrm{BH}}\sim 3.16\times {10}^8{\mathrm{M}}_{\odot }$ [111] and the resulting Eddington luminosity is $L_{\mathrm{Edd}}\sim 3.98\times {10}^{46}$ erg/s. Although jet loading at super-Eddington rates may occur during brief episodes of flaring activity, it is unrealistic to expect such conditions to persist during extended periods of steady, quiescent emission [110]. For both sources the injection of protons follows a simple power-law distribution with ${\gamma }_{\min }=100$, ${\gamma }_{\max }={10}^6$ and ${\alpha }_p=1$ (spectral index).

Figure 3 and Figure 4 display the multiwavelength SED for the blazars 1ES 0414+009 and 1ES 1959+650 along with observational data from several catalogs, including BeppoSAX, VERITAS, Fermi-LAT and others, covering frequencies from radio to gamma rays. The dark blue and red lines correspond to synchrotron emission and IC scattering of primary electrons (best fitting parameters from JetSeT). The black line represents synchrotron emission and IC scattering of protons. The light blue line represents Bethe-Heitler pair production ($p\gamma \to p{e}^{+}{e}^{-}$) from ${10}^{-3}$ eV to ${10}^{15}$ eV (in the case of 1ES 0414+009) and ${10}^{16}$ eV (in the case of 1ES 1959+650). The green line refers to the synchrotron emission and IC scattering of electron-positron pair production ($\gamma \gamma \to e^{+}{e}^{-}$), which is observed in the energy range from ${10}^{-3}$ eV to ${10}^{15}$ eV (in the case of 1ES 0414+009) and ${10}^{17}$ eV (in the case of 1ES 1959+650). The yellow line corresponds to the multi-wavelength emission from proton-photon interaction generating charged pions ($p\gamma \to {\pi }^{\pm }\to {\mu }^{\pm }\to e^{\pm }$) in the range from ${10}^{-1}$ eV to ${10}^{15}$ eV (in the case of 1ES 0414+009) and ${10}^{17}$ eV (in the case of 1ES 1959+650). The purple dotted curve illustrates the gamma-ray emission resulting from the decay of neutral pions produced in proton–photon interactions ($p\gamma \to {\pi }^0\to \gamma \gamma$) observed in the range from ${10}^{11}$ eV to ${10}^{16}$ eV (in the case of 1ES 0414+009) and ${10}^{10}$ eV to ${10}^{17}$ eV (in the case of 1ES 1959+650). The gamma-ray emission from the decay of neutral pions produced in proton–proton interactions ($pp\to {\pi }^0\to \gamma \gamma$) is illustrated by the dark blue dashed line, which is observed in the range ${10}^{11}$ eV to ${10}^{15}$ eV (in the case of 1ES 0414+009) and ${10}^9$ eV to ${10}^{17}$ eV (in the case of 1ES 1959+650).

In Figure 3, the model suggests a possible hadronic contribution within the energy range of ${10}^2$ to ${10}^4$ eV. In this framework, Bethe–Heitler pair production, resulting from the interaction of ultra-high-energy protons with ambient photon fields, leads to the generation of secondary electron–positron pairs. These secondary pairs emit synchrotron radiation as they spiral around magnetic field lines, producing the observed X-ray fluxes in the model. In addition to synchrotron emission, the same Bethe–Heitler secondary pairs also undergo IC scattering, transferring energy to low-energy ambient photons and boosting them into the gamma-ray regime. This IC emission gives rise to the high-energy bump seen in the blue curve, which accurately reproduces the gamma-ray flux of $1.78\times {10}^{-11}$ erg ${\mathrm{cm}}^{-2}$ ${\mathrm{s}}^{-1}$ at $1.08\times {10}^{12}$ eV, as measured by the Fermi-LAT and reported in the 2FHL catalog [88]. Furthermore, the high-energy flux of $3.2\times {10}^{-12}$ erg ${\mathrm{cm}}^{-2}$ ${\mathrm{s}}^{-1}$ at approximately ${10}^{12}$ eV, detected by the ARGO-YBJ observatory (ARGO2LAC catalog) [112], is attributed to IC scattering by secondary electron–positron pairs produced through $\gamma \gamma$ interactions. Thus, both the synchrotron and IC emissions from Bethe–Heitler secondaries, combined with IC scattering from $\gamma \gamma$-induced pairs, contribute significantly to the broadband emission of the source, offering a coherent and self-consistent interpretation of the multiwavelength observations.

In Figure 4, the model indicates a potential hadronic contribution that fills the "gap" between the two characteristic broadband emission features. In this scenario, Bethe–Heitler pair production by ultra-high-energy protons accounts for the X-ray fluxes observed in this intermediate region. These X-ray emissions originate from synchrotron radiation produced by secondary electron–positron pairs generated through the Bethe–Heitler process. Additionally, the IC component, represented by the high-energy bump in the same blue curve, contributes significantly to explaining the observational data between approximately ${10}^{11}$ eV and ${10}^{14}$ eV. This energy range corresponds to VHE gamma rays detected during flaring activity by the Whipple observatory. Notably, on June 4, 2002, the source exhibited a dramatic gamma-ray flare without a simultaneous increase in the X-ray band, marking the first clear detection of an “orphan” gamma-ray flare from a blazar [96]. The analysis of such sources underscores the necessity of adopting a lepto-hadronic framework to fully account for their high-energy emission behavior. This approach has important implications, as it suggests the potential for future neutrino detections and provides strong support for the existence of nuclear acceleration processes in these environments [113,114,115,116].

5. Conclusions

We employed the open-source softwares JetSeT and ${\mathrm{AM}}^3$ to model the SEDs of two high-frequency-peaked BL Lac (HBL) sources: 1ES 0414+009 and 1ES 1959+650. We determined the best-fit model parameters by matching the predicted multi-wavelength emission to publicly available observational data for both sources. Since the overlapping of different instruments, some points have multiple values. So we averaged our multiband data in this work for the fitting process. The derived results have indicated that a purely leptonic model is insufficient to explain the high-energy gamma-ray emissions observed in both sources. Therefore, we used a lepto-hadronic model and analyzed the interactions between high-energy protons and ambient photons within blazar jets to investigate gamma-ray emission from charged particles. The inclusion of hadronic interactions ($pp$ and $p\gamma$) significantly improves the fit of the model for the two HBLs, suggesting a substantial contribution of hadronic processes to their high-energy emission in a flare state. For both sources, the Bethe-Heitler pair production ($p\gamma \to p+e^{+}+e^{-}$) by protons describes the X-ray and very high-energy gamma-ray fluxes. Bethe-Heitler pair production injects secondary electrons and positrons into the radiation zone. These secondaries spiral in magnetic fields (producing synchrotron radiation) and scatter soft photons via IC, contributing to the high-energy part of the SED. This process complements the emission from primary electrons. In conclusion, the lepto-hadronic framework effectively reproduces the multiwavelength observations across a broad energy spectrum for the two sources, underscoring the importance of incorporating hadronic processes to achieve a comprehensive understanding of BL Lac emission.