From Gamma Rays to Cosmic Rays: Lepto-Hadronic Modeling of Blazar Sources as Candidates for Ultra-High-Energy Cosmic Rays

1. Introduction

Cosmic rays (CRs) with energies above ${10}^{19}$ eV have Larmor radii that surpass the scale of the Milky Way, indicating that they likely originate from extragalactic sources. Data from the Pierre Auger Observatory suggest a possible correlation between the arrival directions of the highest-energy cosmic rays (UHECRs, $E>{10}^{18}$ eV) and nearby active galactic nuclei (AGN) [1]. However, this connection remains uncertain and requires further observational support to be confirmed. The foundational idea for particle acceleration in astrophysical environments was introduced by Enrico Fermi in 1949 [2], who proposed that charged particles could gain energy through stochastic interactions with magnetic fields. This mechanism is applicable to acceleration processes in shock fronts, such as those formed in the aftermath of gravitational collapses or around accreting black holes at the centers of galaxies [3], including those found in blazars [4].

Blazars have been extensively observed across multiple wavelengths, particularly in the radio, X-ray, and gamma-ray bands, and are known for their variability in both flux and spectral features on timescales ranging from several years down to mere minutes [5,6,7]. These sources have also been proposed as potential emitters of high-energy astrophysical neutrinos [6,8,9]. With the accumulation of long-term observational data, especially from the IceCube Neutrino Observatory, statistical searches for neutrino counterparts have become increasingly viable [10].

Understanding the mechanisms responsible for the broadband emission observed in blazars is crucial for interpreting their multi-messenger signatures. Their electromagnetic radiation is predominantly non-thermal, and the spectral energy distributions (SEDs) typically exhibit two broad peaks. The low-energy peak, spanning the infrared to X-ray range, is widely attributed to synchrotron radiation from relativistic electrons within the jet. Depending on the peak frequency of this synchrotron component, BL Lac objects are further classified as low-, intermediate-, or high-frequency peaked (LBL, IBL, and HBL, respectively). The high-energy peak, located in the gamma-ray band, remains the subject of ongoing investigation. Leptonic models attribute this emission to inverse Compton scattering, involving either synchrotron photons from the jet itself (in synchrotron self-Compton, SSC, models) or external photon fields such as those from the accretion disk, BLR, or dusty torus (in external Compton, EC, models). Hadronic scenarios, on the other hand, invoke processes such as synchrotron emission from ultra-relativistic protons and photohadronic interactions, which may also lead to the production of high-energy neutrinos [8,9] (see also [11] and references therein).

In this study, we perform a detailed SED analysis of three blazar sources: AP Librae, 1H 1914–194, and PKS 0735+178. These sources are of particular interest due to their detection in the gamma-ray band and, in some cases, at very-high-energy (VHE) gamma rays, offering a valuable opportunity to probe the physical conditions within their relativistic jets. Through broadband multiwavelength modeling, we aim to constrain key jet parameters such as magnetic field strength, emission region size, particle energy distributions, and bulk motion. Additionally, we assess the relative contributions of leptonic and hadronic emission processes to the high-energy component of their SEDs. This approach enables us to investigate the potential role of these blazars as sources of ultra-high-energy cosmic rays and high-energy neutrinos, contributing to the broader understanding of particle acceleration and multi-messenger astrophysics in AGN jets.

The structure of the paper is as follows: in Section 2, we describe the selection of our blazar sample and summarize their key observational properties. Section 3 presents the modeling framework, the computational tools employed, and the methodology used for parameter estimation. This includes the physical assumptions adopted, the spectral fitting procedures, and the criteria for model selection. Section 4 presents the results of our modeling along with a discussion of the implications for blazar physics and multi-messenger astrophysics. We conclude our findings in Section 5. Throughout this work, we adopt a flat $\Lambda$CDM cosmology with the Hubble constant $H_0=67.77\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, matter density parameter ${\Omega }_M=0.27$, radiation density ${\Omega }_r=0$, and dark energy density ${\Omega }_{\Lambda }=0.73$.

2. Blazar Sample Overview: AP Librae, 1H 1914–194, and PKS 0735+178

In this section, we present an overview of the three blazar sources selected for leptonic and lepto-hadronic modeling: AP Librae, 1H 1914–194, and PKS 0735+178. These sources were chosen based on their gamma-ray detection, multiwavelength coverage, and potential relevance as candidates for UHECR production.

2.1. AP Librae

AP Librae, one of the earliest identified BL Lac objects [12], has a reliably measured redshift of $z=0.049\pm 0.002$ [13], consistent with the 6dF galaxy survey measurement [14]. The source was initially detected in the radio band as PKS 1514–24 during the Parkes survey [15], and later associated with the optically variable counterpart AP Librae [16,17]. Its host galaxy contains a supermassive black hole with a mass of approximately ${10}^{8.4}{M}_{\odot }$ based on stellar velocity dispersion measurements [18].

In X-rays, the source was first detected by the Einstein Observatory [19]. At high energies ($E>100$ MeV), it was tentatively identified as 3EG J1517−2538 in the EGRET catalog, though its soft photon index ($\Gamma =2.66\pm 0.43$) implied a low flux at very high energies (VHE) [20]. Earlier VHE observations using the Mark 6 gamma-ray telescope yielded only upper limits [21,22].

AP Librae was later confirmed as a gamma-ray emitter by Fermi-LAT (0FGL J1517.9−2423) with a harder photon index ($\Gamma =1.94\pm 0.14$) [23], increasing the likelihood of detection at VHE. In 2010, H.E.S.S. announced the first VHE detection of the source [24]. Subsequent studies presented a broadband SED from radio to TeV [25], and Chandra observations revealed an extended X-ray jet [26], making AP Librae the only known TeV-detected BL Lac with a resolved X-ray jet.

The high-energy to very-high-energy (HE–VHE) emission of AP Librae has been modeled using a log-parabolic function, yielding a peak energy of approximately 450 MeV (${log}_{10}(E_{\mathrm{peak}}/\mathrm{MeV})=2.65\pm 0.{93}_{\mathrm{stat}}\pm 0.{45}_{\mathrm{sys}}$) [27]. This value lies near the lower bound of the Fermi-LAT sensitivity and may be interpreted as an upper limit. Such a low peak energy is uncommon among extragalactic VHE emitters, which typically exhibit peak emissions above hundreds of GeV. The broadband SED of AP Librae is similarly atypical, being dominated by inverse-Compton emission with an X-ray component that cannot be attributed to synchrotron processes. This suggests a common origin for the X-ray and gamma-ray emissions, with a high-energy component that may be shifted toward lower energies, potentially even below the Fermi-LAT range [25,28].

The VHE emission from AP Librae has been explored through a lepto-hadronic framework in which the extended kiloparsec-scale jet serves as the primary emission region [29]. In this model, synchrotron radiation from ultrarelativistic protons is proposed as the mechanism behind the observed TeV gamma rays. To match the VHE spectrum detected by H.E.S.S., the model requires extreme conditions, including proton energies exceeding $3.98\times {10}^{21}$ eV and magnetic fields around 1 mG. Such parameters lead to a total jet power estimate of $\sim 5\times {10}^{48}$ erg s⁻¹, which significantly exceeds the Eddington luminosity for the source. Despite this, the model successfully reproduces the broadband SED when combined with standard synchrotron and SSC emission from electrons in the inner jet. However, the study concludes that while proton synchrotron emission is a viable explanation, the extreme energetics required impose strong constraints on the plausibility of this scenario [29].

2.2. 1H 1914–194

1H 1914–194, located at a redshift of $z=0.137$ [30], was first detected by NASA’s HEAO 1 (High-Energy Astronomy Observatory 1) during its mission from 1977 to 1979. The primary objective of HEAO 1 was to survey the sky for X-ray and gamma-ray point sources across the 150 eV to 10 MeV energy range, and it was notably the first mission to detect diffuse background radiation (see [31]). Among the sources in the current sample, 1H 1914–194 remains undetected in the TeV gamma-ray band and has the most limited multiwavelength coverage. Detailed radio observations of this X-ray-selected BL Lacertae object (XBL) were later provided by Kollgaard et al. (1996) [31], who described it as having a bright core and a jet-like feature extending northwest, ending in a localized brightening within a diffuse halo.

1H 1914–194 is listed in the Third Catalog of Hard Fermi-LAT Sources (3FHL), which compiles results from the first seven years of Fermi-LAT observations using the Pass 8 data analysis. Within the energy range of 10–1000 GeV, the source exhibits a photon flux of approximately 3.96 × 10⁻¹⁰ photons cm⁻² s⁻¹, with a statistical uncertainty of 3.70 × 10⁻¹¹. Additionally, a power-law fit across the broader 10 MeV to 1000 GeV range yields an energy flux of 2.78 × 10⁻¹¹ erg cm⁻² s⁻¹, accompanied by an uncertainty of 4.45 × 10⁻¹² [32].

1H 1914–194 has garnered attention as a potential candidate in high-energy neutrino studies. The IceCube Neutrino Observatory detected a flux of diffuse cosmic neutrinos with energies in the PeV range during 2013 and 2014 [33], but the sources of these neutrinos remain unidentified. Among the proposed astrophysical counterparts are HBLs, which have been suggested as plausible neutrino emitters due to the potential for hadronic acceleration within their jets to energies exceeding ${10}^{17}$ eV, adequate for neutrino production through hadronic interactions [34,35]. Interestingly, while HBLs show some degree of spatial correlation with the IceCube neutrino arrival directions, other blazar subclasses, such as Flat Spectrum Radio Quasars (FSRQs), despite their higher gamma-ray luminosity and richer external photon environments, appear not to be associated with these events [33,36] . This observation challenges some expectations and further highlights the importance of investigating individual BL Lac sources such as 1H 1914–194.

While Resconi et al. (2017) [37] report a statistical significance of 2.91$\sigma$ (p-value of 0.18%) for a possible association between blazar subclasses and high-energy neutrinos, the current angular resolution of IceCube events does not allow for a definitive identification of individual astrophysical sources. Despite this limitation, 1H 1914–194 has emerged as a promising candidate, with its flaring activity and gamma-ray emission aligning with the arrival directions of some detected neutrinos [34]. Based on this association, the expected contribution from this source to the IceCube signal is estimated at approximately 0.44 muon neutrinos per year [34]. Among the potential counterparts studied, Petropoulou et al. (2015) [33] identify 1H 1914–194 as having the strongest consistency between the modeled neutrino flux and IceCube observations.

2.3. PKS 0735+178

PKS0735+178, classified as an intermediate-synchrotron-peaked BL Lac (IHBL), is among the most luminous BL Lac objects across multiple wavelengths. In the NRAO VLA Sky Survey (NVSS), its radio flux density at 1.4 GHz reaches 2.3 Jy, ranking it as the fifth brightest BL Lac in the 5th edition of the Roma-BZCat catalog [38,39]. Historical radio observations from the early 1990s reported flux densities as high as ∼ 5 Jy at 4.8 and 8 GHz, placing the source among the most powerful blazars overall [40]. In the high-energy (HE; $E>100$ MeV) gamma-ray band, PKS0735+178 also stands out, being listed among the top 20 most luminous IHBLs in the Fermi 4LAC-DR2 catalog [41].

The optical spectrum of the source is notably featureless, complicating the determination of a precise redshift. A strong intervening absorption line has established a firm lower limit of $z\ge 0.424$ [42,43,44], while the tentative detection of its host galaxy led to an estimated redshift of $z=0.45\pm 0.06$ [45]. More recently, a redshift of $z\sim 0.65$ has been proposed based on the hypothesis that the source is associated with a group of nearby faint galaxies [46]. Even adopting the conservative estimate of $z=0.424$, the inferred radio and gamma-ray luminosities are exceptionally high for this class of objects, with $L_{\mathrm{R}}\sim {10}^{27}$ WHz⁻¹ and $L_{\gamma }\sim {10}^{47}$ erg s⁻¹.

Notably, during its intense flaring activity in December 2021, PKS 0735+178 was temporally and spatially correlated with high-energy neutrino detections reported by multiple observatories, including IceCube [47], Baikal-GVD [48], Baksan [49], and KM3NeT [50], highlighting its relevance in multi-messenger astrophysics.

Further support for this potential association was provided by Prince et al. [51], who conducted a detailed lepto-hadronic modeling of the blazar’s broadband emission to explore its role as a neutrino source. They focused on the ∼172 TeV IceCube track-like event, in addition to the Baikal and KM3NeT detections, and tested scenarios with and without external blackbody photon fields as targets for photohadronic interactions. Although the predicted neutrino fluxes fall below IceCube’s current sensitivity, the analysis revealed a strong connection between the X-ray emission and possible neutrino production. Their findings indicate that PKS 0735+178 may contribute non-negligibly to the diffuse astrophysical neutrino background. As emphasized in their study, future neutrino detectors with enhanced angular resolution and sensitivity will be crucial to identifying such sources and refining theoretical flux predictions.

3. Multi-wavelength SED Modeling of BL Lacs

3.1. Leptonic SED fitting

For the leptonic modeling of the SEDs, we used version 1.3.0 of the JetSeT codebase [52,53,54], an open-source platform implemented in C and Python. JetSeT is designed to simulate both the radiation and particle dynamics in relativistic jets, applicable to a wide range of astrophysical objects. It enables detailed modeling of several leptonic emission processes, including synchrotron, synchrotron self-Compton (SSC), and inverse Compton scattering of external photon fields originating from structures such as the accretion disk, broad-line region (BLR), dusty torus (DT), and the cosmic microwave background (CMB). In addition, it accounts for attenuation due to pair production with EBL photons, using widely accepted models [55,56,57].

Our modeling adopts a one-zone leptonic scenario where the emission region is approximated as a homogeneous sphere. The size of this region is linked to the observed variability timescale through the relation

$$R={\displaystyle \frac{c{t}_{\mathrm{var}}\delta }{1+z}},$$

(1)

with $t_{\mathrm{var}}=1\mathrm{day}$ [58,59]. The emitting blob travels with a bulk Lorentz factor $\Gamma$ and is viewed at an angle $\theta$, resulting in a Doppler factor:

$$\delta ={\left[\Gamma (1-\beta cos\theta )\right]}^{-1}$$

(2)

[60]. The electron population within the jet is assumed to follow a broken power-law distribution [61,62,63]:

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

(3)

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

JetSeT’s modeling strategy involves an initial exploratory phase followed by a detailed physical fit. The first step leverages the SEDShape module, which performs empirical fits using power-law and log-parabolic functions across the observed SED. These empirical profiles help to approximate the synchrotron and SSC components, thereby informing the parameter ranges used in the subsequent modeling. Physical modeling is then initialized via the ObsConstrain module, using the broken power-law form in Eq. 3.

The parameter space is explored using a Bayesian framework, where prior distributions are imposed to constrain the models to physically meaningful values. JetSeT expresses all parameters using Astropy units, ensuring interoperability with other Python-based libraries. The fitting routines support both frequentist and Bayesian methodologies. For deterministic optimization, JetSeT integrates with iminuit [64] and SciPy’s bounded least-squares solver [65]. For probabilistic inference, we employed the MCMC engine provided by emcee [66]. Throughout the modeling, the redshift is fixed, and we assume a cold proton content corresponding to 10% of the relativistic electron energy density, consistent with [67].

The code then calculates the kinetic power carried by relativistic electrons in the AGN frame as

$$L_{\mathrm{e}}=\pi R^2{\Gamma }^2{\beta }_{\Gamma }c{m}_{\mathrm{e}}{c}^2{\int }_1^{\infty }d{\gamma }_{\mathrm{e}}{n}_{\mathrm{e}}\left({\gamma }_{\mathrm{e}}\right){\gamma }_e$$

(4)

and compares it with the magnetic power (Poynting flux), given by

$$L_{\mathrm{B}}=\pi R^2{\Gamma }^2{\beta }_{\Gamma }c{\displaystyle \frac{B^2}{8\pi }}.$$

(5)

In addition, the code computes the kinetic luminosity of cold protons, under the assumption of charge neutrality provided by one cold proton per electron (i.e., neglecting electron-positron pairs):

$$L_{\mathrm{p},\mathrm{cold}}=\pi R^2{\Gamma }^2{\beta }_{\Gamma }c{n}_{\mathrm{e}}{m}_{\mathrm{p}}{c}^2.$$

(6)

If electron-positron pairs are also present in the jet, they would reduce the total kinetic luminosity, which is defined as $L_{\mathrm{kin}}\equiv L_{\mathrm{e}}+L_{\mathrm{p},\mathrm{cold}}$.

3.2. Lepto-Hadronic SED Modeling

We investigated hadronic scenarios using the open-source AM³ (Astrophysical Multi-Messenger Modeling) framework [68], which self-consistently evolves coupled populations of photons, electrons, positrons, protons, neutrons, and neutrinos (plus intermediate products) in an isotropic magnetic field. AM³ solves the time-dependent kinetic equations governing non-linear processes, such as electromagnetic cascades and secondary-particle feedback, thereby capturing the full lepto-hadronic interplay within astrophysical sources.

Our emission zone is represented by a spherical “blob” of comoving radius $R^{\prime }$, propagating with bulk Lorentz factor $\Gamma$ along the jet. Primary electrons and protons are injected isotropically: electrons follow a broken power-law energy spectrum, while protons obey a single power-law up to a maximum Lorentz factor ${\gamma }_{\mathrm{p},\max }^{\prime }$. High-energy protons interact with the ambient photon bath according to the prescriptions of Hümmer et al. [69], yielding charged and neutral pions. Subsequent pion decays produce secondary gamma rays, neutrinos, electrons, and positrons, which seed further electromagnetic cascades.

AM³ includes all relevant radiative and particle-interaction processes: synchrotron emission (with self-absorption), inverse-Compton scattering by both electrons and protons, Bethe–Heitler pair production

$$p+\gamma \to p+e^{-}+e^{+},$$

photon–photon pair production and annihilation, and the dynamical evolution of the resulting secondaries. The magnetic field $B^{\prime }$ is treated as turbulent and randomly oriented within the blob, significantly influencing the resulting SED and multi-messenger signatures.

By incorporating the Domínguez et al. (2011) [70] extragalactic background light (EBL) attenuation model into our lepto-hadronic framework, we accurately account for photon absorption along the line of sight to Earth. This allows us to quantify the hadronic component of the blazar’s SED while simultaneously modeling both neutrino and gamma-ray emission from a single compact region. The resulting combined fit offers crucial constraints on the temporal and spectral correlations between high-energy neutrino events and gamma-ray flare episodes in blazars.

4. Results and Discussion

Utilizing the approaches described in Section 3.1 and Section 3.2, we carried out a time-dependent, non-thermal spectral modeling of the blazars AP Librae, 1H 1914–194, and PKS 0735+178 across multiple wavebands. To this end, we compiled archival measurements, spanning from radio frequencies through very-high-energy gamma-rays, sourced from the Firmamento database1. Initially, the broad-band spectral points were fitted via a Markov Chain Monte Carlo algorithm implemented in JetSeT, and the resulting posterior parameter distributions were then adopted as inputs for the subsequent lepto-hadronic simulations using the AM³ framework.

4.1. Synchrotron-Self Compton (SSC) Modeling

Each free parameter in our fitting procedure is constrained by physically motivated limits. We fix the source redshift and set the emission-region distance from the central black hole to $R_{\mathrm{H}}={10}^{17}\mathrm{cm}$ (the JetSeT default), and prescribe allowed intervals for the remaining quantities to accelerate convergence and avoid unphysical solutions. In particular, we bound the emission-region radius R, the magnetic field strength B, the minimum, break, and maximum Lorentz factors ${\gamma }_{min}$, ${\gamma }_{\mathrm{b}}$, ${\gamma }_{max}$, as well as the particle spectral indices p and $p_1$. These ranges are chosen to conform with values reported in the literature [71,72]. The resulting best-fit parameter values are listed in Table 1. The calculated energy densities and luminosities for the three sources are presented in Table 2. The total radiative output is given by

$$L_{\mathrm{rad}}\approx L_{\mathrm{syn}}+L_{\mathrm{SSC}},$$

(7)

where $L_{\mathrm{syn}}$ and $L_{\mathrm{SSC}}$ represent the powers emitted via synchrotron and synchrotron self-Compton processes, respectively [73].

Figure 1, Figure 2, and Figure 3 present the one-zone SSC fits to the SEDs of AP Librae, 1H 1914–194, and PKS 0735+17, respectively. In AP Librae (Figure 1), the synchrotron peak occurs at $\sim 1.5\times {10}^{14}$ Hz, characteristic of an intermediate-synchrotron-peaked (ISP) BL Lac, while the inverse-Compton peak appears at $\sim 1.0\times {10}^{22}$ Hz, extending into the very-high-energy gamma-ray band. In contrast, both 1H 1914–194 (Figure 2) and PKS 0735+178 (Figure 3) exhibit high-synchrotron-peaked (HSP) BL Lac behavior, with synchrotron peaks at $\sim 9.5\times {10}^{15}$ Hz and $\sim 3.0\times {10}^{15}$ Hz, and inverse-Compton peaks at $\sim 1.0\times {10}^{25}$ Hz and $\sim 2.0\times {10}^{24}$ Hz, respectively. In all three cases, the residuals plotted in the lower panels cluster tightly around zero, indicating that the SSC model provides a robust fit to the multiwavelength observations.

4.2. Lepto-Hadronic Spectral Energy Distribution Modeling

The lepto–hadronic SED modeling was performed using the publicly available AM³ code. The nonthermal electron population was adopted from our JetSeT fits (see Section 4.1), yielding jet electron luminosities of

$$L_{\mathrm{e}}\approx 1.03\times {10}^{45}\mathrm{erg}{\mathrm{s}}^{-1}$$

for AP Librae,

$$1.14\times {10}^{45}\mathrm{erg}{\mathrm{s}}^{-1}$$

for 1H 1914–194, and

$$2.65\times {10}^{45}\mathrm{erg}{\mathrm{s}}^{-1}$$

for PKS 0735+17. For the hadronic component, the proton luminosity $L_{\mathrm{p}}$ was constrained to the Eddington limit,

$$L_{\mathrm{p}}\le L_{\mathrm{Edd}},$$

where

$$L_{\mathrm{Edd}}=1.3\times {10}^{38}{\displaystyle \frac{M_{\mathrm{BH}}}{M_{\odot }}}\mathrm{erg}{\mathrm{s}}^{-1}.$$

Adopting black hole masses of $M_{\mathrm{BH}}\approx 2.5\times {10}^8{M}_{\odot }$ for AP Librae [74], $5\times {10}^8{M}_{\odot }$ for 1H 1914–194 [75], and $2.5\times {10}^9{M}_{\odot }$ for PKS 0735+178 [76] yields Eddington luminosities of

$$L_{\mathrm{Edd}}\approx 3.25\times {10}^{46},8.19\times {10}^{46},\mathrm{and}3.15\times {10}^{47}\mathrm{erg}{\mathrm{s}}^{-1},$$

respectively. While transient super–Eddington jet loading may occur during flares, sustained quiescent emission at such rates is unlikely [77]. Informed by a suite of relativistic shock acceleration simulations, which employ both Monte Carlo methods and particle–in–cell techniques to model proton energization at collisionless shocks, spectral indices in the range ${\alpha }_{\mathrm{p}}\sim 1.2\text{–}1.6$ are routinely recovered for strong, ultrarelativistic shocks with magnetic obliquities up to ${60}^{\circ }$ (e.g., [78,79]). These studies demonstrate that the combination of rapid first–order Fermi acceleration and efficient particle scattering off self–generated turbulence leads to unusually hard power–law tails compared to the canonical $p\approx 2$ result for nonrelativistic shocks. Motivated by these findings, we adopt a single power–law injection spectrum for the jet protons of the form

$${\displaystyle \frac{dN}{d\gamma }}\propto {\gamma }^{-{\alpha }_{\mathrm{p}}},$$

with

$${\gamma }_{\min }=100,{\gamma }_{\max }={10}^6,{\alpha }_{\mathrm{p}}=1.5.$$

Here, ${\gamma }_{\min }=100$ represents the approximate threshold above which protons begin to participate in the diffusive acceleration cycle, while ${\gamma }_{\max }={10}^6$ is set by balancing the acceleration timescale against both radiative losses and particle escape in the jet environment.

Rodrigues et al. (2024) [77] adopt a conservative strategy in which the X-ray and radio measurements serve only as upper bounds on the leptonic contribution. However, for some objects, a standard electron synchrotron model can in fact account for the observed flux in these bands. This choice stems from the compactness of the high-energy emission zone, which produces both gamma-rays and neutrinos, and its concomitant synchrotron self-absorption at radio frequencies, rendering it optically thick below GHz wavelengths [77]. Extending this hadronic picture, we incorporate proton–proton ($pp$) collisions alongside the originally considered photohadronic ($p\gamma$) interactions. In $pp$ reactions, a relativistic proton impacts a thermal proton, producing both neutral and charged pions that subsequently decay via

$${\pi }^0\to \gamma \gamma ,{\pi }^{+}\to {\mu }^{+}+{\nu }_{\mu },{\mu }^{+}\to e^{+}+{\nu }_e+{\overline{\nu }}_{\mu },$$

thus providing additional channels for very-high-energy gamma-ray and neutrino generation. By unifying $p\gamma$ and $pp$ processes, including all relevant flavor and charge channels, our model offers a more complete depiction of hadronic interactions and their role in shaping the multiwavelength and neutrino outputs of blazar jets.

Figure 4, Figure 5, and Figure 6 display the multiwavelength SEDs of AP Librae, 1H 1914–194, and PKS 0735+17, respectively, together with archival data from radio through gamma-rays. The red curve corresponds to the synchrotron and inverse-Compton emission of primary electrons (JetSeT best-fit); the black curve shows proton synchrotron and proton-induced inverse-Compton. The light blue curve denotes Bethe–Heitler pair production ($p\gamma \to p{e}^{+}{e}^{-}$), while the green curve marks the emission from secondary $e^{+}{e}^{-}$ pairs produced via photon–photon interactions ($\gamma \gamma \to e^{+}{e}^{-}$). The yellow curve represents the cascade from charged-pion decay ($p\gamma \to {\pi }^{\pm }\to {\mu }^{\pm }\to e^{\pm }$), and the purple dotted line illustrates gamma-rays from neutral-pion decay in photopion processes ($p\gamma \to {\pi }^0\to \gamma \gamma$). The dark blue dashed curve shows the gamma-ray output from neutral-pion decay in proton–proton collisions ($pp\to {\pi }^0\to \gamma \gamma$). Finally, the pink curve represents the all-flavor neutrino fluxes from $p\gamma$ interactions, while the brown curve from $pp$ interactions.

In Figure 4 and Figure 6, our lepto–hadronic model naturally produces an emission component that bridges the gap between the low-energy synchrotron peak and the high-energy inverse-Compton peak: ultra-high-energy protons interact with ambient photons via the Bethe-Heitler process ($p+\gamma \to p+e^{+}+e^{-}$), injecting secondary electrons and positrons whose synchrotron radiation reproduces the observed X-ray flux in the intermediate band ($\sim {10}^3$–${10}^4$ eV). Simultaneously, charged pions from photopion interactions ($p+\gamma \to {\pi }^{\pm }$) decay through ${\pi }^{\pm }\to {\mu }^{\pm }\to e^{\pm }$, and these decay electrons upscatter target photons via inverse-Compton to form the pronounced “yellow” hump in the ${10}^{10}$–${10}^{11}$ eV range, thereby filling the SED valley and achieving excellent agreement with the multiwavelength data. Furthermore, as shown in Figure 5, the same population of electrons, both those injected by Bethe–Heitler pair production and those originating from muon decay, can inverse-Compton scatter soft photons up to very high energies (${10}^{10}$–${10}^{12}$ eV and beyond), accounting for the observed VHE gamma-ray tail in that source.

5. Conclusions

In this work, we have carried out a comprehensive multiwavelength study of three blazar sources, AP Librae, 1H 1914–194, and PKS 0735+178, using both leptonic (one-zone SSC) and time-dependent lepto-hadronic (AM³) models. Our SSC fits reproduce the low- and high-energy peaks of the SEDs with characteristic synchrotron and inverse-Compton components, confirming AP Librae as an ISP BL Lac and 1H 1914–194 and PKS 0735+178 as HSP BL Lacs. By coupling these results with lepto-hadronic simulations, we demonstrate that secondary processes, Bethe–Heitler pair cascades and photopion interactions, naturally fill the SED “valley,” accounting for the observed X-ray emission, and extend the gamma-ray spectrum into the VHE band.

We find that proton injection with ${\alpha }_p=1.5$, ${\gamma }_{min}=100$, and ${\gamma }_{max}={10}^6$ remains energetically viable under $L_p\le L_{\mathrm{Edd}}$, while providing adequate neutrino fluxes in both $p\gamma$ and $pp$ channels. Notably, 1H 1914–194 exhibits an inverse–Compton tail from secondary electrons that can account for photons up to ${10}^{12}$ eV. Our modeling yields total jet kinetic powers of $\sim {10}^{46}$–${10}^{47}$ erg s⁻¹ and electron–to–proton energy ratios consistent with other blazar studies.

These findings support the scenario in which HSP BL Lacs serve as promising UHECR accelerators and neutrino emitters. Future observations with next-generation Cherenkov telescopes and high-resolution neutrino detectors will be critical for testing the hadronic signatures predicted here.