EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 710 (June 2026) A&A, 710 (2026) A279 Full HTML
Open Access
Issue
A&A
Volume 710, June 2026
Article Number A279
Number of page(s) 8
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202558081
Published online 19 June 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

Relativistic jets launched by active galactic nuclei (AGNs) emit radiation across the entire electromagnetic spectrum, from radio waves to gamma rays, and possibly also produce neutrinos. These jets provide unique laboratories for studying black hole physics, particle acceleration, and high-energy emission processes (e.g., Blandford et al. 2019). Blazars are among the best targets for investigating jet physics. They are radio-loud AGNs whose relativistic jet is oriented close to our line of sight. Because of this alignment, the nonthermal jet emission is strongly amplified by relativistic beaming and dominates the spectral energy distribution (SED) (e.g., Romero et al. 2017; Böttcher 2019).

The SED of a blazar typically exhibits two broad humps. The low-energy hump peaks between the infrared and X-ray bands and is generally attributed to synchrotron emission from relativistic leptons. The high-energy hump, which peaks from MeV to TeV energies, is less well understood, and its origin remains debated. In the leptonic scenario, it is explained as the result of inverse-Compton scattering, where relativistic leptons upscatter photons from either their own synchrotron radiation, also called synchrotron self-Compton (SSC), or from external radiation fields (e.g., Maraschi et al. 1992; Sikora et al. 1994; Ghisellini et al. 1998). Alternatively, hadronic models attribute the high-energy component to processes involving relativistic protons, such as direct proton synchrotron emission or synchrotron radiation from secondary particles produced in hadronic interactions (e.g., Cerruti 2020; Sol & Zech 2022). The hadronic picture is supported by the possible association of some blazars with high-energy neutrino detections (e.g., Aartsen 2018).

The positions of the two SED peaks are used to classify blazars, giving rise to the so-called blazar sequence, in which the peak frequencies are anticorrelated with the bolometric luminosity (Fossati et al. 1998; Donato et al. 2001; Ghisellini et al. 2017). At the low-luminosity end of this sequence lie extreme high-frequency-peaked BL Lacs (EHBLs), which are the most efficient particle accelerators among blazars (for a review, Biteau et al. 2020).

In the seminal paper by Ghisellini (1999), the author proposed extending the blazar sequence beyond its previously assumed limits, suggesting the existence of a new BL Lac population whose synchrotron spectrum peaks in the MeV band. In principle, there is no strict theoretical bound on the maximum synchrotron photon energy, except for the so-called synchrotron burn-off limit. If nonthermal particles are accelerated at shocks (or, more generally, by a gyroresonant process), the minimum acceleration timescale corresponds to the gyration time, tg = 2πrg/c, where rg is the particle gyroradius. This process is constrained by synchrotron cooling, which imposes an upper limit on the achievable photon energy that is independent of the magnetic field strength. For a relativistically beamed source, the maximum synchrotron photon energy can reach ∼25 δ MeV, comfortably exceeding the energies observed in EHBLs (e.g., Guilbert et al. 1983; de Jager et al. 1996). Specifically for high-energy BL Lacs, assuming that the acceleration timescale is a multiple of the gyration time, tacc = ηtg, leads to very high values, η ∼ 105, implying that shocks are rather slow accelerators (e.g., Inoue & Takahara 1996; Garson et al. 2010).

A natural starting point for modeling MeV-peaking BL Lacs is one of the competing frameworks proposed for EHBLs. Owing to their peculiar spectral features, namely the large separation between the two peaks and their steep gamma-ray spectrum, the standard one-zone leptonic model typically requires unusually low magnetic fields (B < 10 mG) and very high average electron Lorentz factors ( γ ¯ ∼ 10 3 − 10 4 Mathematical equation: $ \bar{\gamma} \sim 10^3{-}10^4 $). These values place EHBLs as clear outliers with respect to the parameter ranges inferred for other BL Lacs (Tavecchio et al. 2010; Kaufmann et al. 2011; Costamante et al. 2018). To account for these differences, several alternative scenarios have been proposed, including a Maxwellian-like electron distribution (Lefa et al. 2011), a beam of high-energy hadrons (Essey & Kusenko 2010), emission from a large-scale jet (Aharonian et al. 2008), lepto-hadronic models (Cerruti et al. 2015), and multiple shock acceleration (Zech & Lemoine 2021).

In Sciaccaluga & Tavecchio (2022) and Sciaccaluga et al. (2024), we presented a novel framework for EHBLs based on the combined action of shock and turbulence acceleration. In this scenario, nonthermal leptons are initially accelerated by a shock and are subsequently energized by the downstream turbulence.

All EHBL models can, in principle, be extended to even higher energies, potentially giving rise to a new population of BL Lacs. If these extreme sources are indeed observed, they would provide a unique laboratory for testing competing acceleration and emission scenarios.

BL Lacs whose synchrotron hump peaks in the MeV band, here termed ultra-extreme high-frequency-peaked BL Lacs (UEHBLs), would be ideal targets for the upcoming missions designed to explore this energy range (for a review about extragalactic sources, see Sbarrato et al. 2025). Historically, the 0.1–100 MeV band has been poorly sampled: only a few missions have partially covered it and did so with limited sensitivity. For this reason, this portion of the electromagnetic spectrum is commonly referred to as the MeV gap. In recent years, several satellites have been proposed to bridge this gap, spanning small-, medium-, and large-scale missions, such as the Compton Spectrometer and Imager (COSI, Tomsick 2022), All-sky Medium Energy Gamma-ray Observatory eXplorer (AMEGO-X, Caputo et al. 2022), and enhanced ASTROGAM (e-ASTROGAM, de Angelis et al. 2018).

If UEHBLs do exist, they may have already left observable signatures in neighboring energy bands. For instance, in the hard X-rays, satellites such as Swift Burst Alert Telescope (BAT, Barthelmy et al. 2005) provide excellent sky coverage and sensitivity. Interestingly, the BAT catalogs include a number of sources without clear counterparts, which might represent promising candidates for this hypothesized BL Lac population (Oh et al. 2018; Lien et al. 2025).

The paper is organized as follows. In Section 2 we describe the hybrid shock–turbulence acceleration model. Section 3 reports and discusses three representative realizations of the model. In Section 4 we examine the potential temporal and polarimetric characteristics of UEHBLs and discuss the implications for possible acceleration and emission mechanisms. Throughout the paper, the following cosmological parameters are assumed: H0 = 70 km s−1 Mpc−1, ΩM = 0.3, and ΩΛ = 0.7.

2. Shock-turbulence model

In this section, we briefly describe our shock-turbulence acceleration model. This is one possible EHBL model that might be extended to UEHBLs (for more details, see Sciaccaluga & Tavecchio 2022 and Sciaccaluga et al. 2024). As mentioned in the introduction, Zech & Lemoine (2021) proposed a model based on multiple shock acceleration. Their idea relied on the fact that recollimation (or, more generally, standing) shocks occur in series, as demonstrated by several 2D fluid simulations (e.g., Gomez et al. 1995; Mizuno et al. 2015). However, 3D fluid simulations revealed that the flow is subject to instabilities that evolve into turbulence and ultimately disrupt the cycle of recollimation and reflection shocks (e.g., Matsumoto & Masada 2013; Gourgouliatos & Komissarov 2018; Boula et al. 2025; Hu et al. 2025; Costa et al. 2026). For this reason, we assumed that particles are first accelerated by a shock and then further energized by the downstream turbulence.

We assumed that the emitting region is the turbulent downstream of a recollimation shock, and we adopted a leaky-box spatially averaged approach, so that our scenario effectively was a one-zone leptonic model. The temporal evolution of the electrons and turbulence is governed by two coupled Fokker–Planck equations (e.g., Eilek 1979; Miller et al. 1996; Kakuwa 2016; Gong et al. 2025),

{ ∂ f ∂ t = 1 p 2 ∂ ∂ p [ p 2 D p ∂ f ∂ p + p 2 ( ∂ p ∂ t ) rad f ] − f t esc + I f ∂ W ∂ t = ∂ ∂ k [ k 2 D k ∂ ∂ k ( W k 2 ) ] − W t dam + I w , Mathematical equation: $$ \begin{aligned} \left\{ \begin{aligned}&\frac{\partial f}{\partial t} = \frac{1}{p^2} \frac{\partial }{\partial p} \left[ p^2 D_p \frac{\partial f}{\partial p} + p^2 \left( \frac{\partial p}{\partial t} \right)_{{\mathrm{rad}}} f \right] - \frac{f}{t_{\mathrm{esc}}} +I_f \\&\frac{\partial W}{\partial t} = \frac{\partial }{\partial k} \left[ k^2 D_k \frac{\partial }{\partial k} \left(\frac{W}{k^2} \right) \right] - \frac{W}{t_{\mathrm{dam}}} + I_w \end{aligned}\right., \end{aligned} $$(1)

where p is the electron momentum, k the fluctuation wavenumber, f(p, t) the electron isotropic phase-space density, and W(k, t) the turbulence energy density per unit of wavenumber. Eq. (1) describes the interplay between particle processes (resonant acceleration, cooling, escape, and injection) and turbulence processes (cascading, damping, and injection). All quantities are defined in the comoving frame of the emission region. The expressions for particle acceleration and turbulence damping are derived from quasi-linear theory, although alternative approaches exist (e.g., Lemoine et al. 2024).

The diffusion coefficient of electrons is given by

D p = p 2 β a 2 c U B r g 2 ∫ k res W B k d k , Mathematical equation: $$ \begin{aligned} D_p = \frac{p^2 \,\beta _a^2 \, c}{U_B \, r_g^2} \int _{k_{\mathrm{res}}} \frac{W_B}{k}\,dk, \end{aligned} $$(2)

where βa is the dimensionless Alfvén speed, UB = B2/8π is the magnetic energy density, rg is the electron gyroradius, WB ≈ W/2 is the magnetic component of the turbulence energy spectrum, and kres = 1/rg is the resonant wavenumber.

The electron cooling term accounts for synchrotron and inverse-Compton contributions, following standard formulae from the literature (e.g., Jones 1965; Chiaberge & Ghisellini 1999).

The electron escape time is equal to

t esc = R c + R 2 κ ‖ , Mathematical equation: $$ \begin{aligned} t_{\mathrm{esc}} = \frac{R}{c} + \frac{R^2}{\kappa _\parallel }, \end{aligned} $$(3)

where R is the emission region radius, and κ∥ = crg/9ζ(kres) is the spatial diffusion coefficient along the magnetic field, with ζ(k) = kWB/UB the relative amplitude of the turbulent magnetic field energy density for a given k. When the mean free path is large, the escape time is essentially the geometric escape time. Conversely, when the mean free path is small, turbulence traps particles in the emission region, leading to a longer escape time.

Before gaining energy through turbulence, particles are initially accelerated at the shock, which thus acts as the injector for the emission region. The electron injection number density per unit of Lorentz factor In is given by

I n = I n , 0 γ − p e − γ γ cut with γ > γ min , Mathematical equation: $$ \begin{aligned} I_n = I_{n,0} \,\gamma ^{-p}\, e^{-\frac{\gamma }{\gamma _{\mathrm{cut}}}} \quad {\mathrm{with}}\ \, \gamma > \gamma _\mathrm{min} , \end{aligned} $$(4)

where γ denotes the electron Lorentz factor, In, 0 is the injection normalization, p = 2 is the power-law slope, and γmin = 103 and γcut = 105 are the minimum and cutoff Lorentz factor, respectively, of the injected electron distribution. We note that In = 4πp2mec If, where If is the corresponding injection distribution in phase space. The slope and cutoff of the injection were determined from recent simulations of diffusive shock acceleration for weakly magnetized shocks (e.g., Vanthieghem et al. 2020; Zech & Lemoine 2021). The injection was normalized to the injected electron power Pn,

P n = V ∫ γ m e c 2 I n d γ , Mathematical equation: $$ \begin{aligned} P_n = V \int \gamma \, m_e\, c^2\, I_n \,d\gamma , \end{aligned} $$(5)

where V = 10πR3 is the emission region volume, modeled as a cylinder with radius R and length 10 R. The length of the emission volume, related to the region where the instability develops and triggers turbulence in the plasma, was roughly estimated on fluid simulations (e.g., Matsumoto et al. 2021; Boula et al. 2025; Costa et al. 2026). We expect that within this distance, the magnetic field decay and the adiabatic losses effectively quench the emission.

In a Kolmogorov phenomenology, the diffusion coefficient of turbulence is given by

D k = k 3 β a c kW 2 U B · Mathematical equation: $$ \begin{aligned} D_k = k^3\, \beta _a\, c \sqrt{\frac{k W}{2 U_B}}\cdot \end{aligned} $$(6)

Without strong damping, and given this diffusion coefficient with continuous injection, W(k) would evolve toward the standard Kolmogorov spectrum, W(k)∝k−5/3 (Zhou & Matthaeus 1990).

The turbulence damping time is equal to

t dam = − 4 π e 2 β a 2 m e c k ∫ γ res γ 2 ∂ ∂ γ ( n e γ 2 ) d γ , Mathematical equation: $$ \begin{aligned} t_{\mathrm{dam}} = -\frac{4\pi e^2 \beta _a^2 }{m_e c k} \int _{\gamma _{\mathrm{res}}} \gamma ^2 \frac{\partial }{\partial \gamma } \left( \frac{n_e}{\gamma ^2} \right) d\gamma , \end{aligned} $$(7)

where ne(γ) = 4πmecp2f(p) is the electron number density per unit of Lorentz factor, and γres is the resonant Lorentz factor, defined as k = 1/rg(γres). Damping is determined by enforcing energy conservation, so that the energy driving electron acceleration is removed from the turbulence.

The turbulence injection term is equal to

I W = I W , 0 δ ( k − k 0 ) , Mathematical equation: $$ \begin{aligned} I_W = I_{W,0}\, \delta (k-k_0), \end{aligned} $$(8)

where δ is the Dirac function, Iw, 0 is the normalization, and k0 = 1/L is the injection wavenumber, with L = R/10. The injection is normalized to the injected turbulence power Pw,

P w = V ∫ I W d k . Mathematical equation: $$ \begin{aligned} P_w = V \int I_W \,dk. \end{aligned} $$(9)

The radiative output of the emission region was computed using the SSC model, following the standard formulae in the literature (e.g., Jones 1968; Blumenthal & Gould 1970; Ghisellini et al. 1988). For the treatment of relativistic beaming, we adopted the usual blob amplification formula (see Appendix C of Zech & Lemoine 2021).

We employed the Chang–Cooper algorithm (Chang & Cooper 1970) on a logarithmic grid, using 20 points per decade for momentum, wavenumber, and frequency. The system was evolved over a time interval of 10 R/c with 100 time steps.

3. Results

In this section, we discuss some examples of realizations of the hybrid shock–turbulence model that produce BL Lacs with a spectral peak in the MeV band. In Fig. 1, we present three representative realizations. The model depends on six parameters: the emission region radius R, the dimensionless Alfvén speed βa, the magnetic field strength B, the injected electron power Pe, the injected turbulence power Pw, and the relativistic Doppler factor 𝒟 = [Γ(1 − β cos θv)]−1 (where Γ and β denote the bulk Lorentz factor and the dimensionless velocity of the fluid, respectively, while θv is the observer’s viewing angle). The parameter sets corresponding to the three realizations are summarized in Table 1. In all cases, the redshift was fixed at z = 0.14, matching that of the prototypical EHBL 1ES 0229+200 (Woo et al. 2005). In addition to the input parameters, the table also lists two derived quantities we used to check the model consistency: the magnetization, σ = βa2/(1 − βa2), and the relative amplitude of turbulent magnetic fluctuations, δB/B = (∫WB dk/UB)1/2.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Spectral energy densities for three model realizations: Case A (dashed red), case B (dash–dotted blue), and case C (dash–double–dotted green), including attenuation from extragalactic background light absorption (Franceschini & Rodighiero 2017). The solid dark and light red lines show COSI 2- and 9-year sensitivities (Tomsick 2022), and the dark yellow and orange lines show AMEGO-X (3 yr; Caputo et al. 2022) and e-ASTROGAM (1 yr; de Angelis et al. 2018) sensitivities. The Fermi-LAT 10-year extragalactic sensitivity (taken from https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm) is shown as a solid blue line, and the CTA North and South 50-hour sensitivities (taken from https://www.ctao.org/for-scientists/performance/) are shown as dark and light green lines. The red shaded region marks the BAT-band flux of Swift J1949.7–3636 (Lien et al. 2025), a candidate MeV-peaking BL Lac. For reference, an SED template of a giant elliptical galaxy taken from Silva et al. (1998) is displayed in dotted grey, renormalized to the magnitude of the host galaxy of 1ES 0229+200 (Costamante et al. 2018).

Table 1.

Physical parameters and derived quantities for the three modeled cases.

In all three cases, the emission region radius was fixed at subparsec scales, as expected for BL Lacs (Tavecchio et al. 1998, 2010). The dimensionless Alfvén speed was set to ensure low magnetization, which makes shock acceleration efficient and favors the development of instabilities and turbulence in the downstream region (e.g., Vanthieghem et al. 2020; Matsumoto et al. 2021). The remaining parameters were varied as needed. However, since we adopted the quasi-linear approximation, we verified that the relative turbulence amplitude remained small, tha tis, δB/B ≪ 1. The time evolution of the electron and turbulence spectra for case A is reported in Appendix A. Cases B and C exhibit a similar behavior and are therefore not shown.

Case A provides the most efficient acceleration, with the low-energy bump peaking at ∼2 MeV. In the MeV band, the low-energy peak is fully detectable by AMEGO-X and e-ASTROGAM, with the latter also able to detect case A beyond the peak. COSI cannot observe the source within its nominal mission lifetime, namely 2 years, although case A would become detectable if the mission were extended. Detection by COSI would also be possible in the event of a strong flare. Although EHBLs, the closest known class of sources to those hypothesized here, do not exhibit strong flaring at the frequencies corresponding to the low-energy bump, radiative thermal instabilities could trigger such flares, as discussed in Section 4.

The relative turbulence amplitude in case A is close to the nominal boundary of the quasilinear regime. However, it is possible to obtain spectra with the same synchrotron peak energy for lower levels of magnetic turbulence, that is, lower values of δB/B. As discussed in Sciaccaluga et al. (2024), the model parameter space is significantly degenerate, which allows different combinations of physical quantities to yield similar spectral outcomes. For example, we can increase the Alfvén speed parameter βa while reducing the injected turbulence power (see case D in Appendix C). This leads to a higher magnetization, with corresponding values as high as σ ∼ 0.1, and a reduced turbulence level δB/B ∼ 0.1. However, this relatively high magnetization might suppress the efficiency of shock acceleration and limit the development of downstream turbulence. Alternatively, a decrease in the injected turbulence power can be compensated for by an increase in the Doppler factor. This allows for configurations with δB/B ∼ 0.1 that still produce a spectrum peaking in the same region as case A (see case E in Appendix C). These examples illustrate that the appearance of spectra is not uniquely tied to a narrow region of parameter space, but instead arises from a broader set of physically plausible configurations. A comprehensive and systematic exploration of the parameter degeneracy is beyond the scope of the present work and will be addressed in future studies.

Potential constraints might arise from the hard X-ray band. As mentioned in the introduction, valuable information is expected from Swift-BAT, which operates in the hard X-rays and offers excellent sky coverage. Recently published BAT catalogs (Oh et al. 2018; Lien et al. 2025) contain several candidates that may represent UEHBLs. To identify such candidates, we filtered the BAT catalogs with the following criteria: (i) integrated flux in the BAT band above 10−11 erg cm−2 s−1, (ii) photon index Γ < 1.5, (iii) likely extragalactic origin (defined by a Galactic latitude |b|> 10°), (iv) absence of a clear counterpart in the soft X-ray and optical bands. The requirements on the flux and index in the hard X-rays were imposed to identify sources that might be detectable by the next generation of MeV telescopes. Applying these requirements, we identified ten possible sources according to the catalog classifications. Their number might drop to six when possible X-ray counterpart candidates are updated (see Appendix B for further details). In Fig. 1, we plot the integrated BAT flux of one of the candidates, SWIFT J1949.7-3636. All BAT candidates lack a soft X-ray counterpart that would allow us to clearly associate them with sources in optical catalogs. This is consistent with case A, in which the soft X-ray flux is close to the sensitivity limit of current soft X-ray satellites. A potential approach to identifying a soft X-ray counterpart would be to perform a deep observation of the BAT error region using a soft X-ray telescope with a sufficiently large field of view, such as Swift X-Ray Telescope (XRT, Burrows et al. 2005) and XMM-Newton (Jansen et al. 2001). In the event of multiple detections, spectral characteristics (e.g., a hard photon index) might help us to distinguish among the possible candidates. These candidates could be FSRQs with their high-energy peak located in the MeV band. In this scenario, BAT would be detecting the rising part of the high-energy component. However, UEHBLs and MeV FSRQs are easily distinguishable. First, given a soft X-ray detection providing better localization, the optical counterpart should present broad emission lines in the case of FSRQs. Second, MeV FSRQs after the peak typically exhibit a power-law tail extending into the GeV band, which should be detectable by Fermi Large Area Telescope (LAT).

The case A parameters were selected to ensure consistency with the BAT data, but this requirement is highly restrictive. For realizations B and C, we therefore chose not to impose this constraint. From a theoretical perspective, there is no strict limitation on the properties of potential MeV BL Lacs, and intermediate cases between EHBLs and case A should, in principle, also be possible. Moreover, as previously highlighted, these sources may not remain stable over time. Case A is consistent with a BAT source that has been monitored for 157 months. However, during this period, the source may have experienced flaring episodes, during which its brightness exceeded the average flux measured by BAT.

Case B has a low-energy bump peaking at ∼0.5 MeV, and it is the brightest candidate in the MeV band, potentially detectable by COSI within its nominal mission lifetime. Case C, with a peak at 0.2 MeV, represents the scenario most closely resembling EHBLs. While COSI could still detect case C within its lifetime, the signal would lie near the edge of its sensitivity range, allowing observation only just beyond the peak. The superior sensitivities of AMEGO-X and e-ASTROGAM would enable them to probe the peak and higher energies in cases B and C.

None of the three cases is observable by Fermi or the Cherenkov Telescope Array (CTA) because the scattering cross section in the Klein–Nishina regime is severely suppressed in a full leptonic scenario (see Section 4 for further details on a possible hadronic component).

Finally, all three cases, the optical band is dominated by the galaxy emission, even more than in EHBLs, as shown in Fig. 1. A more accurate localization, for example, using soft X-ray observations, would enable us to detect the host galaxy of these sources.

4. Discussion

We demonstrated that the hybrid shock–turbulence scenario, one of the possible models for EHBLs, can also account for UEHBLs. We presented three representative model realizations and discussed their detectability with proposed MeV observatories and with existing telescopes at other wavelengths. In addition, we identified ten potential candidates from the BAT catalogs.

Proposed MeV satellites such as COSI, AMEGO-X, and e-ASTROGAM are essential for detecting these new sources. Their observations would provide valuable insights into particle acceleration. Previous studies have already shown that MeV observations of powerful FSRQ (sources at the opposite extreme of the blazar sequence with respect to EHBLs) could reveal the presence of a thermal component in the particle spectrum in addition to the standard nonthermal power law. This would indicate that shocks are the underlying acceleration mechanism (Tavecchio et al. 2025). On the other hand, the detection of BL Lacs peaking in the MeV band would imply that electrons can achieve extremely high Lorentz factors. This condition might help us to further distinguish among the different particle acceleration mechanisms proposed in the literature. Conversely, a non-detection of such sources could be interpreted in two ways. First, particle acceleration, regardless of the underlying mechanism, must proceed relatively slowly, as discussed in Section 1. Second, UEHBLs may simply be too faint to be detected by the proposed MeV missions.

In addition to the spectral properties of UEHBLs, we are also interested in their polarimetric and temporal features. In the MeV band, we expect relatively high polarization degrees (≳40%), following the trend observed in HBLs and EHBLs (e.g., Di Gesu et al. 2022; Liodakis et al. 2022; Ehlert et al. 2023; Kouch et al. 2024), where the polarization is strongly chromatic, that is, the degree of polarization increases with frequency. Based on the flux predicted by our model (excluding strong flaring states), the polarization of UEHBLs would remain undetectable for COSI and AMEGO-X, while e-ASTROGAM might measure it, depending on the actual polarization degree.

It is interesting to note that UEHBLs might represent an exceptional class of sources. In contrast to standard blazars, their emission would be dominated by synchrotron radiation, since inverse-Compton contribution is strongly quenched by Klein-Nishina suppression. This peculiarity might affect the variability properties of these sources. A source dominated by the pressure of nonthermal electrons and whose cooling is dominated by synchrotron losses might indeed be subject to the radiative thermal instability discussed by Marscher (1980). A small increase in the magnetic field in limited portions of the flow can, in fact, lead to a collapse of the region due to the increased cooling rate of the relativistic electrons and consequent loss in pressure. The rapid radiative losses suffered by the electrons during the compression would produce a rapid increase in the emissivity, with flares potentially detectable in the MeV band.

As already outlined, in a purely leptonic scenario, the inverse-Compton emission is expected to be suppressed, since most electrons scatter photons in the Klein–Nishina regime. Consequently, as illustrated in Fig. 1, UEHBLs are expected to remain undetectable by CTA. However, we cannot exclude that other mechanisms can result in a detectable flux even at these energies. First, the emission might originate from a nonthermal hadronic component producing synchrotron radiation at TeV energies (e.g., Mannheim 1993; Aharonian 2000; Mücke & Protheroe 2001). Alternatively, a second emission region farther downstream in the jet might exist, where nonthermal electrons generate infrared/radio photons. These photons might then be upscattered by the nonthermal electrons in the primary emission region in the Thomson regime, giving rise to a bright inverse-Compton bump detectable by CTA. Variability studies might help us to distinguish between these two scenarios. A rapid variability on timescales of hours would disfavor the hadronic model, since the synchrotron cooling timescale for protons is expected to be orders of magnitude longer.

Acknowledgments

This work has been funded by ASI under contract 2024-11-HH.0. We acknowledge financial support from an INAF Theory Grant 2024 (PI F. Tavecchio) and the European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (2022C9TNNX).

References

  1. Aharonian, F. A. 2000, New Astron., 5, 377 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aharonian, F. A., Khangulyan, D., & Costamante, L. 2008, MNRAS, 387, 1206 [NASA ADS] [CrossRef] [Google Scholar]
  3. Barthelmy, S. D., Barbier, L. M., Cummings, J. R., et al. 2005, Space Sci. Rev., 120, 143 [Google Scholar]
  4. Biteau, J., Prandini, E., Costamante, L., et al. 2020, Nat. Astron., 4, 124 [Google Scholar]
  5. Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
  6. Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [Google Scholar]
  7. Böttcher, M. 2019, Galaxies, 7, 20 [Google Scholar]
  8. Boula, S., Tavecchio, F., Bodo, G., et al. 2025, A&A, 704, A200 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165 [Google Scholar]
  10. Caputo, R., Ajello, M., Kierans, C. A., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 044003 [NASA ADS] [CrossRef] [Google Scholar]
  11. Cerruti, M. 2020, Galaxies, 8, 72 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910 [Google Scholar]
  13. Chang, J. S., & Cooper, G. 1970, J. Comput. Phys., 6, 1 [Google Scholar]
  14. Chiaberge, M., & Ghisellini, G. 1999, MNRAS, 306, 551 [Google Scholar]
  15. Costa, A., Bodo, G., Tavecchio, F., et al. 2026, A&A, 705, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Costamante, L., Bonnoli, G., Tavecchio, F., et al. 2018, MNRAS, 477, 4257 [NASA ADS] [CrossRef] [Google Scholar]
  17. de Angelis, A., Tatischeff, V., Grenier, I. A., et al. 2018, J. High Energy Astrophys., 19, 1 [NASA ADS] [CrossRef] [Google Scholar]
  18. de Jager, O. C., Harding, A. K., Michelson, P. F., et al. 1996, ApJ, 457, 253 [Google Scholar]
  19. Di Gesu, L., Donnarumma, I., Tavecchio, F., et al. 2022, ApJ, 938, L7 [CrossRef] [Google Scholar]
  20. Donato, D., Ghisellini, G., Tagliaferri, G., & Fossati, G. 2001, A&A, 375, 739 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Ehlert, S. R., Liodakis, I., Middei, R., et al. 2023, ApJ, 959, 61 [NASA ADS] [CrossRef] [Google Scholar]
  22. Eilek, J. A. 1979, ApJ, 230, 373 [CrossRef] [Google Scholar]
  23. Essey, W., & Kusenko, A. 2010, Astropart. Phys., 33, 81 [Google Scholar]
  24. Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433 [Google Scholar]
  25. Franceschini, A., & Rodighiero, G. 2017, A&A, 603, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Garson, A. B., III, Baring, M. G., & Krawczynski, H. 2010, ApJ, 722, 358 [Google Scholar]
  27. Ghisellini, G. 1999, Astropart. Phys., 11, 11 [NASA ADS] [CrossRef] [Google Scholar]
  28. Ghisellini, G., Guilbert, P. W., & Svensson, R. 1988, ApJ, 334, L5 [NASA ADS] [CrossRef] [Google Scholar]
  29. Ghisellini, G., Celotti, A., Fossati, G., Maraschi, L., & Comastri, A. 1998, MNRAS, 301, 451 [NASA ADS] [CrossRef] [Google Scholar]
  30. Ghisellini, G., Righi, C., Costamante, L., & Tavecchio, F. 2017, MNRAS, 469, 255 [NASA ADS] [CrossRef] [Google Scholar]
  31. Gomez, J. L., Marti, J. M. A., Marscher, A. P., Ibanez, J. M. A., & Marcaide, J. M. 1995, ApJ, 449, L19 [Google Scholar]
  32. Gong, X.-W., Liu, R.-Y., Zhang, Z.-L., Asano, K., & Lemoine, M. 2025, ApJ, 989, 99 [Google Scholar]
  33. Gourgouliatos, K. N., & Komissarov, S. S. 2018, Nat. Astron., 2, 167 [Google Scholar]
  34. Guilbert, P. W., Fabian, A. C., & Rees, M. J. 1983, MNRAS, 205, 593 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hu, X.-F., Mizuno, Y., & Fromm, C. M. 2025, A&A, 693, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. IceCube Collaboration (Aartsen, M. G., et al.) 2018, Science, 361, eaat1378 [NASA ADS] [Google Scholar]
  37. Inoue, S., & Takahara, F. 1996, ApJ, 463, 555 [Google Scholar]
  38. Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Jones, F. C. 1965, Phys. Rev., 137, 1306 [Google Scholar]
  40. Jones, F. C. 1968, Phys. Rev., 167, 1159 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kakuwa, J. 2016, ApJ, 816, 24 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kaufmann, S., Wagner, S. J., Tibolla, O., & Hauser, M. 2011, A&A, 534, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Kouch, P. M., Liodakis, I., Middei, R., et al. 2024, A&A, 689, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Lefa, E., Rieger, F. M., & Aharonian, F. 2011, ApJ, 740, 64 [NASA ADS] [CrossRef] [Google Scholar]
  45. Lemoine, M., Murase, K., & Rieger, F. 2024, Phys. Rev. D, 109, 063006 [Google Scholar]
  46. Lien, A. Y., Krimm, H. A., Markwardt, C. B., et al. 2025, ApJ, 989, 161 [Google Scholar]
  47. Liodakis, I., Marscher, A. P., Agudo, I., et al. 2022, Nature, 611, 677 [CrossRef] [Google Scholar]
  48. Mannheim, K. 1993, A&A, 269, 67 [NASA ADS] [Google Scholar]
  49. Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5 [CrossRef] [Google Scholar]
  50. Marscher, A. P. 1980, ApJ, 239, 296 [Google Scholar]
  51. Matsumoto, J., & Masada, Y. 2013, ApJ, 772, L1 [NASA ADS] [CrossRef] [Google Scholar]
  52. Matsumoto, J., Komissarov, S. S., & Gourgouliatos, K. N. 2021, MNRAS, 503, 4918 [NASA ADS] [CrossRef] [Google Scholar]
  53. Miller, J. A., Larosa, T. N., & Moore, R. L. 1996, ApJ, 461, 445 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38 [Google Scholar]
  55. Mücke, A., & Protheroe, R. J. 2001, Astropart. Phys., 15, 121 [Google Scholar]
  56. Oh, K., Koss, M., Markwardt, C. B., et al. 2018, ApJS, 235, 4 [Google Scholar]
  57. Romero, G. E., Boettcher, M., Markoff, S., & Tavecchio, F. 2017, Space Sci. Rev., 207, 5 [NASA ADS] [CrossRef] [Google Scholar]
  58. Sbarrato, T., Ajello, M., Buson, S., et al. 2025, Space Sci. Rev., 221, 62 [Google Scholar]
  59. Sciaccaluga, A., & Tavecchio, F. 2022, MNRAS, 517, 2502 [NASA ADS] [CrossRef] [Google Scholar]
  60. Sciaccaluga, A., Tavecchio, F., Landoni, M., & Costa, A. 2024, A&A, 687, A247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 [Google Scholar]
  62. Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103 [Google Scholar]
  63. Sol, H., & Zech, A. 2022, Galaxies, 10, 105 [NASA ADS] [CrossRef] [Google Scholar]
  64. Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 [Google Scholar]
  65. Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570 [Google Scholar]
  66. Tavecchio, F., Nava, L., Sciaccaluga, A., & Coppi, P. 2025, A&A, 694, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Tomsick, J.& COSI Collaboration 2022, in 37th International Cosmic Ray Conference, 652 [Google Scholar]
  68. Vanthieghem, A., Lemoine, M., Plotnikov, I., et al. 2020, Galaxies, 8, 33 [NASA ADS] [CrossRef] [Google Scholar]
  69. Woo, J.-H., Urry, C. M., van der Marel, R. P., Lira, P., & Maza, J. 2005, ApJ, 631, 762 [Google Scholar]
  70. Zech, A., & Lemoine, M. 2021, A&A, 654, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Zhou, Y., & Matthaeus, W. H. 1990, J. Geophys. Res., 95, 14881 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Electron and turbulence spectra

Figures A.1 and A.2 show the time evolution of the electron and turbulence spectra of Case A, respectively. Electrons, injected into the downstream region with Lorentz factors in the range 103 < γ ≲ 105, are accelerated by the turbulence up to γ ≳ 107, where radiative cooling becomes dominant. During the acceleration, electrons extract energy from the turbulence, leading to strong damping at wavenumbers corresponding to the Lorentz factors where acceleration is most efficient. At lower wavenumbers, where the turbulence cascade dominates over electron damping, turbulence follows the standard Kolmogorov spectrum. The damping has further consequences. Initially, electrons tend to accumulate at high Lorentz factors. However, as the turbulence is progressively damped, the acceleration efficiency decreases. Consequently, at later times, the injected population becomes dominant again, flattening the spectrum at low Lorentz factors.

Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Time evolution of the electron number density per unit of Lorentz factor as a function of Lorentz factor. The top axis reports the corresponding resonant wavenumber.

Thumbnail: Fig. A.2. Refer to the following caption and surrounding text. Fig. A.2.

Time evolution of the turbulence energy density per unit of wavenumber as a function of wavenumber. The top axis reports the corresponding resonant Lorentz factor.

Appendix B: BAT candidates

We further checked the existence of soft X-ray counterparts for these 10 sources using archival images and catalogs provided by Swift-XRT, Chandra, XMM-Newton, ROSAT, and eROSITA. We confirm the absence of any likely X-ray counterpart for all sources except 4, which show a non-negligible association probability with sources characterized by X-ray fluxes close to the limiting value defined by Lien et al. (2025). Specifically, SWIFT J0007.8-4133 lies 8 arcmin from MCG-07-01-011, a Seyfert 2 galaxy detected by XMM-Newton, Swift-XRT, and eROSITA; SWIFT J0656.0-6560 is 6 arcmin from Fairall 0265, a Seyfert 1 galaxy detected by ROSAT and Swift-XRT; SWIFT J0243.2-0553 is located 2.3 arcmin from the known γ-ray emitting blazar PKS 0240-060; and SWIFT J1949.7-3636 is close to a ROSAT-detected source with no identified optical counterpart.

Table B.1.

Candidate BAT sources with their corresponding right ascension (RA), declination (Dec), observed flux, photon index, and catalogs association class.

Appendix C: Additional cases

Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

The parameters of Case D are R = 1.0 × 1016 cm, B = 1.7 × 10−2 G, βa = 3.0 × 10−1, Pn = 1.0 × 1037 erg/s, Pw = 6.7 × 1039 erg/s, 𝒟 = 1.6 × 101. The magnetization is σ = 1.0 × 10−1, while the relative amplitude of turbulent magnetic fluctuations is δB/B = 1.5 × 10−1

Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

The parameters of Case E are R = 1.0 × 1016 cm, B = 1.7 × 10−2 G, βa = 2.0 × 10−1, Pn = 1.0 × 1037 erg/s, Pw = 2.9 × 1039 erg/s, 𝒟 = 2.0 × 101. The magnetization is σ = 4.2 × 10−2, while the relative amplitude of turbulent magnetic fluctuations is δB/B = 1.3 × 10−1

All Tables

Table 1.

Physical parameters and derived quantities for the three modeled cases.

In the text
Table B.1.

Candidate BAT sources with their corresponding right ascension (RA), declination (Dec), observed flux, photon index, and catalogs association class.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Spectral energy densities for three model realizations: Case A (dashed red), case B (dash–dotted blue), and case C (dash–double–dotted green), including attenuation from extragalactic background light absorption (Franceschini & Rodighiero 2017). The solid dark and light red lines show COSI 2- and 9-year sensitivities (Tomsick 2022), and the dark yellow and orange lines show AMEGO-X (3 yr; Caputo et al. 2022) and e-ASTROGAM (1 yr; de Angelis et al. 2018) sensitivities. The Fermi-LAT 10-year extragalactic sensitivity (taken from https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm) is shown as a solid blue line, and the CTA North and South 50-hour sensitivities (taken from https://www.ctao.org/for-scientists/performance/) are shown as dark and light green lines. The red shaded region marks the BAT-band flux of Swift J1949.7–3636 (Lien et al. 2025), a candidate MeV-peaking BL Lac. For reference, an SED template of a giant elliptical galaxy taken from Silva et al. (1998) is displayed in dotted grey, renormalized to the magnitude of the host galaxy of 1ES 0229+200 (Costamante et al. 2018).
In the text
Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Time evolution of the electron number density per unit of Lorentz factor as a function of Lorentz factor. The top axis reports the corresponding resonant wavenumber.
In the text
Thumbnail: Fig. A.2. Refer to the following caption and surrounding text. Fig. A.2.

Time evolution of the turbulence energy density per unit of wavenumber as a function of wavenumber. The top axis reports the corresponding resonant Lorentz factor.
In the text
Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

The parameters of Case D are R = 1.0 × 1016 cm, B = 1.7 × 10−2 G, βa = 3.0 × 10−1, Pn = 1.0 × 1037 erg/s, Pw = 6.7 × 1039 erg/s, 𝒟 = 1.6 × 101. The magnetization is σ = 1.0 × 10−1, while the relative amplitude of turbulent magnetic fluctuations is δB/B = 1.5 × 10−1
In the text
Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

The parameters of Case E are R = 1.0 × 1016 cm, B = 1.7 × 10−2 G, βa = 2.0 × 10−1, Pn = 1.0 × 1037 erg/s, Pw = 2.9 × 1039 erg/s, 𝒟 = 2.0 × 101. The magnetization is σ = 4.2 × 10−2, while the relative amplitude of turbulent magnetic fluctuations is δB/B = 1.3 × 10−1
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.