© The Authors 2024

1. Introduction

Active galactic nuclei (AGNs) are the most powerful persistent sources in the Universe. The output of these objects is powered by the gravitational energy released by gas accreting onto a supermassive black hole residing in a galactic nucleus. Radio-loud AGNs are characterized by the presence of a relativistic jet produced in the supermassive black hole vicinity that can propagate up to hundreds of kiloparsecs in the most powerful sources. Radio-loud AGNs are further divided into several classes based on their observational features, which are strongly dependent on the observational angle (Urry & Padovani 1995). Blazars are radio-loud AGNs whose jet is pointing toward the observer, and therefore the non-thermal emission of the jet dominates the emission thanks to relativistic Doppler beaming (Romero et al. 2017; Blandford et al. 2019).

The spectral energy distribution (SED) of blazars presents two broad humps, the first due to synchrotron emission by non-thermal electrons. The origin of the second hump, peaking in the gamma-ray band, is still disputed: it could be generated by the interaction of non-thermal electrons with the synchrotron photons or with photons filling the external environment (i.e., synchrotron self-Compton and external Compton models; e.g., Ghisellini et al. 1998) or by hadronic processes, such as proton synchrotron emission or photo-pion production (e.g., Böttcher et al. 2013).

Blazars can be classified using the frequencies of the two peaks (e.g., Ghisellini et al. 2017): the first ranges from the infrared to the X-ray bands and the second from MeV to TeV energies. The most efficient accelerators are the so-called extremely high-frequency-peaked BL Lacs (EHBLs; Costamante et al. 2001), which, in turn, can be divided into two subclasses, extreme synchrotron BL Lacs (which present the first peak above 1 keV) and extreme TeV BL Lacs. In addition, the latter are characterized by the second peak surpassing 1 TeV and a hard GeV spectrum (Γ < 2 with F_ν/ν ∝ ν^−Γ). Furthermore, extreme TeV BL Lacs, at odds with the general behavior of blazars, present a low temporal variability at high energies (Biteau et al. 2020). The SEDs of these sources are difficult to explain through a leptonic single-zone model with a single shock acceleration: the spectral features imply a large minimum Lorentz factor, a small magnetic field far from equipartition, and a power law index incompatible with the theory of diffusive shock acceleration (p < 2 with dN/dE ∝ E^−p). Several solutions have been proposed, such as a Maxwellian-like electron distribution (Lefa et al. 2011), a beam of high-energy hadrons (Essey & Kusenko 2010), internal absorption (Aharonian et al. 2008), emission from a large-scale jet (Böttcher et al. 2008), lepto-hadronic models (Cerruti et al. 2015), and multiple shock acceleration (Zech & Lemoine 2021).

To explain the phenomenology of extreme TeV BL Lacs, we elaborated a double-step acceleration model based on a scenario involving a jet recollimated by the pressure of external material. If the magnetization of the jet is low, which is likely for extreme TeV BL Lacs, the downstream of the shock formed as a consequence of the recollimation becomes unstable and turbulent (Gourgouliatos & Komissarov 2018; Costa et al. 2024). In this scheme, non-thermal particles are first accelerated by the shock and then further energized in the downstream by turbulence through resonant interaction. Comparing turbulence cascading and damping timescales, it is evident that the damping cannot be neglected (Tavecchio et al. 2022). In Sciaccaluga & Tavecchio (2022) we presented a time-dependent one-zone model that includes damping by the accelerating electrons. We calculated the nonthermal electron and turbulence spectra and then derived the SED using the synchrotron self-Compton (SSC) model. We compared our model with data of the prototypical extreme TeV BL Lacs 1ES 0229+200 and adjusted, via visual inspection, our model on the flux points.

In this paper we apply the model sketched above but including other EHBLs with well-sampled SEDs. Moreover, to automatize and parallelize the comparison between the model and data, we developed a procedure based on a Markov chain Monte Carlo (MCMC) sampler. This technique allowed us to explore the parameter space of the model.

The paper is organized as follows: in Sect. 2 we describe the code updates, in Sect. 3 we explain how we use MCMC sampling in our framework, and Sect. 4 is dedicated to the discussion. Throughout the paper, the following cosmological parameters are assumed: H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7.

2. The model

Supposing a proton-electron plasma, the non-thermal electrons are the solely responsible for the emission. As described above, we assumed a model in which particles are accelerated at a shock and then energized through stochastic acceleration in the downstream, where they also emit through synchrotron and SSC mechanism. In our model, the non-thermal population of electrons accelerated at the shock is treated as an injection term in a kinetic equation for the stochastic acceleration.

The time evolution of non-thermal electrons and turbulence is described by a system of two coupled nonlinear Fokker-Planck equations (Eilek 1979; Miller & Roberts 1995; Kakuwa 2016):

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\partial f}{\partial t}=\frac{1}{p^2}\frac{\partial }{\partial p}[p^2{D}_{\mathrm{p}}\frac{\partial f}{\partial p}+p^2{\left(\frac{\partial p}{\partial t}\right)}_{\mathrm{rad}}f]+\frac{f}{t_{\mathrm{esc}}}+I_{\mathrm{f}}\hfill \\ \hfill & \frac{\partial Z}{\partial t}=\frac{1}{k^2}\frac{\partial }{\partial k}\left(k^2{D}_k\frac{\partial Z}{\partial k}\right)+\frac{Z}{t_{\mathrm{dam}}}+\frac{I_{\mathrm{W}}}{k^2}\hfill \end{array}.\end{array}$$(1)

Here f(p, t) is the momentum distribution of non-thermal electrons, and Z(k, t) = W/k², where W(k, t) is the wavenumber energy spectrum of turbulence, that is, the energy density (including both the kinetic and the magnetic component) of waves with wavenumber between k and k + dk. Equation (1) includes all the physical processes related to non-thermal electrons and turbulence: particle resonant acceleration, cooling, escape, and injection together with turbulence cascading, damping, and injection. The two coupled equations are solved numerically using standard schemes (details in Appendix C). We used 20 points per decade for the momentum and wavenumber grid, 50 time steps, and 10 points per decade for the frequency grid.

The momentum diffusion coefficient of the electrons is obtained from the quasi-linear theory of particle-wave interaction, supposing only parallel and antiparallel propagating magnetohydrodynamics waves for a further simplification (e.g., Kakuwa 2016):

$$\begin{array}{c}\hfill D_{\mathrm{p}}=\frac{p^2{v}_{\mathrm{a}}^2}{U_B{r}_{\mathrm{g}}^{2}c}{\displaystyle {\int }_{k_{\mathrm{res}}}}\frac{W_B}{k}\mathrm{d}k,\end{array}$$(2)

where v_a is the Alfvén velocity, U_B = B²/8π the total (ordered plus turbulent) energy density of the magnetic field, m_e the electron mass, c the light speed, r_g = γm_ec/eB the Larmor radius, W_B ≈ W/2 the magnetic field component of the turbulence energy spectrum, and k_res = 2π/r_g is the resonant wavenumber.

In addition to synchrotron cooling, this time we included Compton cooling for electrons, which increased the execution time but, after some optimizations, we were able to reach speeds (approximately a few seconds) comparable with other leptonic codes (e.g., Stathopoulos et al. 2024). For the synchrotron and Compton cooling, we used standard formulae.

In our scenario, the escape depends on the turbulence. When the turbulence is strong, the electrons diffuse in the acceleration region, while if damping is important, the escape time is simply equal to the geometrical escape time. Therefore, the energy-dependent escape timescale can then be written as

$$\begin{array}{c}\hfill t_{\mathrm{esc}}=\frac{R}{c}+\frac{R^2}{{\kappa }_{\Vert }},\end{array}$$(3)

where κ_∥ = cr_g/9ζ(k_res) is the spatial diffusion coefficient along the total magnetic field, while ζ(k) = kW_B/U_B is the relative amplitude of the turbulent magnetic field energy density for a given k.

The injection term describes the particles accelerated at the recollimation shock and advected downstream, supposing a strong shock:

$$\begin{array}{c}\hfill I_n=I_{n,0}{\gamma }^{-2}{e}^{-\frac{\gamma }{{\gamma }_{\mathrm{cut}}}}\mathrm{with}{10}^3<\gamma <{10}^7,\end{array}$$(4)

where I_n is the injected electron number density per unit of time, which can be converted to the injection in the momentum space I_n = 4πp²m_ecI_f. The range of injection is determined from recent simulations of diffusive shock acceleration (Zech & Lemoine 2021). The injection is normalized to the injected non-thermal electron power, P_n:

$$\begin{array}{c}\hfill P_n=V{\displaystyle \int \gamma }{m}_{\mathrm{e}}{c}^2{I}_n\mathrm{d}\gamma ,\end{array}$$(5)

where V = 10πR³ is the emission region volume, modeled as a cylinder with radius R and length 10 R. This estimate of the length for the emission volume, related to the region where the instability develops and triggers turbulence in the plasma, is roughly based on the simulations shown in Matsumoto et al. (2021). We expect that within such a distance, both the magnetic field decay and the adiabatic losses effectively quench the emission (Tavecchio et al. 2022).

Regarding the turbulence equation, we used the Kolmogorov phenomenology, for which the effective diffusion coefficient is

$$\begin{array}{c}\hfill D_k=\frac{1}{2}{k}^3{v}_{\mathrm{a}}\sqrt{\frac{kW}{U_B}}.\end{array}$$(6)

Without strong damping, with this diffusion coefficient and constant injection, W(k) would reach a standard Kolmogorov spectrum, W(k)∝k^−5/3 (Zhou & Matthaeus 1990).

The damping time is obtained by imposing energy conservation, that is, the energy used to accelerate the electrons is subtracted from the turbulence:

$$\begin{array}{c}\hfill t_{\mathrm{dam}}={\left|\frac{4\pi e^2{v}_{\mathrm{a}}^2}{m_{\mathrm{e}}{c}^{3}k}{{\displaystyle \int }}_{\gamma >{\gamma }_{\mathrm{res}}}{\gamma }^2\frac{\partial }{\partial \gamma }\left(\frac{n_{\mathrm{e}}}{{\gamma }^2}\right)\mathrm{d}\gamma \right|}^{-1}.\end{array}$$(7)

Here n_e = 4πm_ecp²f is the electron energy spectrum and γ_res is the resonant Lorentz factor, defined by k = 2π/r_g(γ_res).

Finally, the turbulence is injected at a length L equal to one-tenth of the emitting region radius, from which it cascades to shorter scales:

$$\begin{array}{c}\hfill I_{\mathrm{W}}=I_{\mathrm{W},0}\delta (k-k_0).\end{array}$$(8)

Here δ is the Dirac function and k₀ = 2π/L is the injection wavenumber. The injection is normalized to the injected turbulence power, P_w:

$$\begin{array}{c}\hfill P_{\mathrm{w}}=V{\displaystyle \int I_{\mathrm{W}}}\mathrm{d}k.\end{array}$$(9)

More information on the different terms can be found in Sciaccaluga & Tavecchio (2022).

3. Markov chain Monte Carlo (MCMC)

In the previous exploratory paper (Sciaccaluga & Tavecchio 2022) the model was adjusted on data by eye, a standard procedure in literature. However, this method presents several shortcomings, such as confirmation bias and repetition for each source. Therefore, we decided to move to MCMC sampling. This technique has several advantages: it automatizes and parallelizes the exploration of the parameter space, it can calculate the distribution and the cross-correlation of the model parameters and it permits the imposition of non-diagonal priors (i.e., conditions that depend on multiple parameters).

To explore the parameter space, we used the Python library emcee as the MCMC ensemble sampler (Foreman-Mackey et al. 2013). It implements several moves; we tested the default strech move (Goodman & Weare 2010) and the differential-independence mixture ensemble move (Boehl 2022): with an equal number of steps, the latter presents a lower autocorrelation time, and we therefore decided to adopt it. We assumed a Gaussian likelihood,

$$\begin{array}{c}\hfill ln\mathcal{L}=\frac{1}{2}{\displaystyle \sum _i}\frac{{[{(\nu F)}_{i,m}-{(\nu F)}_i]}^2}{{\sigma }_i^2},\end{array}$$(10)

where (ν F)_i and σ_i are respectively the ith flux point and its uncertainty, while (ν F)_i, m is the model flux at the ith frequency. The model flux is calculated using standard synchrotron and SSC emissivities, after integrating Eq. (1) until t_max = 10R/c (see Sciaccaluga & Tavecchio 2022 for the details). In literature, it is a standard procedure to add a term in the likelihood to account for the scatter due to the non-simultaneity of flux points (Hogg et al. 2010; Stathopoulos et al. 2024), but extreme TeV BL Lacs are highly stable. Tests that included this nuisance parameter resulted in extremely low values; therefore, we neglected it.

3.1. Fixing priors

Our model presents several parameters: the radius of the emission region, R, the Alfvń velocity, v_a, the total magnetic field, B, the non-thermal electron and turbulence injection power in the jet frame, respectively P_n and P_w, and the relativistic Doppler factor, δ.

After some tests, we decided to fix the relativistic Doppler factor to δ = 20. If left free, it tends to reach δ ∼ 100, a value that we consider implausible for a blazar. On the other hand, the relativistic Doppler factor cannot be too small, since it would result in small boosting and in a large emission region radius, again implausible for a blazar region (at sub-parsec scales). We were finally left with five free parameters.

As shown in Table 1, we applied log-flat priors to all parameters over broad ranges, allowing the walkers to explore a large part of the parameter space. Furthermore, we imposed two non-diagonal priors. Since the interaction between turbulence and non-thermal particles is described using quasi-linear theory, we required the turbulent component to be small compared to the total field in the emission region. Moreover, we required the non-thermal component to be a minor part of the post-shock plasma (e.g., Sironi & Spitkovsky 2011). Specifically, we put constraints on the ratio of the energy density of the magnetic turbulence and of the total field, U_B, and on the ratio of the number density of the non-thermal electrons, N_e, and of the thermal plasma, N_p, calculated from the Alfvén velocity and the total magnetic field:

$$\begin{array}{cc}& \frac{\delta B^2}{B^2}=\frac{\int W_B(k)\mathrm{d}k}{U_B}<0.1\hfill \\ \hfill & \frac{N_{\mathrm{e}}}{N_{\mathrm{p}}}<0.1.\hfill \end{array}$$(11)

No strong limits were imposed on the magnetization σ = B²/8πm_pN_p, where m_p is the proton mass. It was used as a consistency check since we expect a low value, necessary for the development of the instabilities in the downstream region.

4. Results

We applied the procedure described above to four extreme TeV BL Lacs taken from Costamante et al. (2018): 1ES 0229+200, 1ES 0347−121, RGB J0710+591, and 1ES 1101−232. All these EHBLs have a well-characterized SED with a synchrotron peak well tracked by Swift and NuSTAR and a high-energy peak covered by Fermi and Cherenkov telescopes. From the sample of Costamante et al. (2018), comprising a total of six sources, we did not consider two EHBLs, namely 1ES 0414+009 and 1ES 1218+304. For the former object, the slope of the Fermi and the TeV spectra are apparently in disagreement, making the SED unsuitable for the automated fit approach. For 1ES 1218+304 the situation is rather complex: in the optical-UV region, the relative contribution of the host galaxy and the jet is not easy to disentangle. Moreover, the Swift data do not clearly trace the shape of the low-energy part of the synchrotron peak. For these reasons, we preferred to exclude also this complex source from our study.

In the X-ray band, the model uncertainty is smaller than in other bands since Swift and NuSTAR measurements are more accurate than Fermi and Cherenkov arrays. It is also worth noticing that the Fermi spectrum is well reproduced for each source, thanks to the softening due to the turbulence damping (Sciaccaluga & Tavecchio 2022). Finally, very high-energy data are compatible within 1σ uncertainty, except for the most energetic flux points of 1ES 0229+200 and 1ES 0347−121, which are still compatible at < 2σ. In no cases did we consider the optical-UV data in the fit, since in these bands the emission is dominated by the host galaxy.

For each source, we initialized 30 walkers, each moving for 10 000 steps and with ∼1000 burn-in steps. In Figs. 1–4 the measured flux data points are shown together with the median and the 1σ credible interval of the calculated SEDs, obtained drawing ∼1000 random samples from the posterior. The corresponding corner plots, showing the distribution of the derived model parameters, are shown in Appendix A.

Fig. 1. Flux points with their errors (black) of 1ES 0229+200 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.

Fig. 2. Flux points with their errors (black) of 1ES 0347−121 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.

Fig. 3. Flux points with their errors (black) of 1ES 1101−232 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.

Fig. 4. Flux points with their errors (black) of RGB J0710+591 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.

Examining the corner plots (see Figs. A.1, A.2, A.3, and A.4), the distributions in logarithmic space appear narrow (evidence of a good reconstruction of the physical parameters), but strongly not-Gaussian. They peak in the same range of parameters, with slightly broader uncertainties for RGB J0710+591 and with 1ES 1101−232 as the only exception. The emitting region radius (R ∼ 10^16 ÷ 17 cm) and the magnetic field (B ∼ 10^{−3 ÷ −2} G) are compatible with the values obtained by phenomenological models that leave the acceleration process unspecified (e.g., Costamante et al. 2018). The Alfvén velocity (v_a ∼ 10⁹ cm s⁻¹) aligns with values obtained in literature (e.g., Kakuwa 2016; Tavecchio et al. 2022). Finally, the injected electron power (P_n ∼ 10^40 ÷ 41 erg s⁻¹) is consistent with values utilized with other BL Lacs (e.g., Ghisellini et al. 2010).

For 1ES 1101−232 the fit provides a radius much larger than the value derived for the other sources (R ∼ 1 pc), together with a smaller magnetic field and higher electron and turbulence powers. These differences are probably connected to the steeply decreasing X-ray spectrum (i.e., a relatively low synchrotron peak frequency), a feature absent in the other sources. A sub-parsec-sized emission region can be obtained by applying a larger relativistic Doppler factor (δ ∼ 40). The results are shown in Appendix B.

Quite interestingly, the corner plots reveal strong cross-correlations among some of the parameters, which are hard to discover through the standard approach based on the visual comparison of the model and data. For example, there is a strong anticorrelation between the radius of the emitting region and the magnetic field for each source. This is explained by the fact that the magnetic field determines the magnetic field energy density U_B, while, for a fixed synchrotron luminosity, the radius regulates the radiation energy density, U_rad. Since the ratio U_B/U_rad is fixed by the ratio of the synchrotron and SSC peak luminosities (e.g., Tavecchio et al. 1998), this fixes the product BR. While this correlation is strong for each source and can be directly interpreted, some other correlations are weaker and much more difficult to explain in simple terms, because they are associated with observational features that involve several parameters.

We calculated the final electron energy distribution, n(γ), and the final turbulence spectrum, W(k), (medians and credible intervals), drafted from the posterior. As a representative case, we display in Figs. 5 and 6 the results for 1ES 0229+200. Finally, we report the median and the credible interval of the final timescales associated with electrons (i.e., acceleration, cooling, and escape time) and turbulence (cascading and damping times) for 1ES 0229+200, respectively in Figs. 7 and 8.

Fig. 5. Final electron energy distribution for 1ES 0229+200 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior.

Fig. 6. Final spectrum of the turbulence for 1ES 0229+200 with 1σ credible intervals (violet) and the median (blue) obtained from the model posterior.

Fig. 7. Final electron timescales (acceleration, cooling, escape) for 1ES 0229+200 with 1σ credible intervals and the median obtained from the model posterior.

Fig. 8. Final turbulence timescales (cascading and damping) for 1ES 0229+200 with 1σ credible intervals and the median obtained from the model posterior.

As discussed in Sciaccaluga & Tavecchio (2022), the combined effects of acceleration and radiative cooling result in a final electron energy distribution displaying a well-defined peak that, for EHBLs, is generally located at Lorentz factors γ_p ≈ 10⁶, where t_acc ∼ t_cool (see Fig. 7). The role of turbulence acceleration is relevant only at high energies, 10⁵ ≲ γ ≲ 10⁶, where it piles up the injected electrons. At a lower energies (γ ≲ 10⁵), particles are inefficiently accelerated and effectively escape from the acceleration and emission region. Since at these energies electrons are unaffected by acceleration or cooling (t_esc ≪ t_cool and t_esc ≪ t_acc), the distribution tends to be similar to that injected by the shock (i.e., γ⁻²).

For turbulence, damping is effective only at large k, where the damping time is less than the cascading time (Figs. 6 and 8). At these wavenumbers, the damping is a consequence of the large number of resonating low-energy electrons. High-energy electrons (γ ≳ 10⁵), responsible for the X-ray emission, resonate with low k modes (k < 10⁻¹¹ cm⁻¹), where, instead, damping is not efficient due to their small number. Notably, the very inefficient acceleration of the electrons at low Lorentz factors at late times (clearly visible in Fig. 7) is related to the strong damping of the turbulence at large k.

The electron energy distribution shows a large uncertainty (about an order of magnitude) linked to the complex interplay between the electron density, the peak energy of the electron energy distribution, the magnetic field and the radius of the emitting region. The uncertainty on the normalization of the electron density, n, can be understood recalling that, as discussed above, the product BR is fixed by the observed ratio of the synchrotron and SSC peak luminosity. This, for a fixed synchrotron luminosity, proportional to L_S ∝ nB²R³, implies that n and R are anticorrelated: nR = const. Therefore, a large (small) radius implies both a low (high) density and a low (high) magnetic field. In order to have a fixed synchrotron peak frequency, proportional to ${\nu }_{\text{s}}\propto {\gamma }_{\text{p}}^{2}B$, the electron distribution must peak at high (low) Lorentz factor γ_p (as visible in Fig. 5). The large uncertainty of the electron distribution is also reflected by the SSC spectrum, which, however, presents a smaller variance. The same SSC spectrum can be obtained by using different electron spectra and target photon fields, therefore reducing the SSC spectrum uncertainty at the expense of less constrained optical and radio spectra.

The turbulence energy spectrum displays a large uncertainty as well. In fact, the same acceleration efficiency can be obtained with a high magnetic field (or Alfvén velocity) and low turbulence injection and vice versa (see Eq. (2) and the corner plot). The much larger uncertainty in the region of k where the spectrum is affected by damping reflects the strong interaction with the electrons.

Starting from the Alfvén velocity, we calculated the magnetization distributions to confirm the consistency of our model. Using 1ES 0229+200 as a representative case, we obtain a low magnetization σ < 10⁻² (see Fig. 9). Such low values are compatible with estimates from simple leptonic scenarios (e.g., Zech & Lemoine 2021) and are below the threshold derived for the development of instabilities in the downstream region of the recollimation shock (Gourgouliatos & Komissarov 2018).

Fig. 9. Distribution of the magnetization of 1ES 0229+200 obtained from the model posterior.

5. Discussion

We have presented the modeling of a sample of EHBLs with well-reconstructed SEDs using a double acceleration (shock+stochastic) model (Tavecchio et al. 2022). By means of the MCMC technique, we have explored the full parameter space. The MCMC sampling permitted us to approach the modeling without confirmation biases and showed us that the different parameters are correlated in a nontrivial way. This confirms that the standard simple visual comparison between the data and one realization of the model is strongly limited and that sampling is necessary to completely characterize the model complexity. Indeed, while the most probable values of the derived model parameters are in line with those derived with the standard approach, in the present work we were able to infer the related uncertainties and the correlations among them. We confirm for all sources that the interplay between the turbulence and the accelerating electrons leads to a strong damping of the turbulent energy at large wave vectors. From the point of view of the energetic budget, the scenario requires moderate powers in the form of injected turbulence and electrons, which is completely compatible with the jet power derived by means of standard one-zone models (Ghisellini et al. 2010).

In the original version of our scenario (Tavecchio et al. 2022) we assumed that the stochastic acceleration occurs in the turbulent region after a recollimation shock. However, constraining the size of the emission region to sub-parsec scales results in a relatively large relativistic factor (δ = 20 − 40). Such values, although modest compared to those derived with simple one-zone leptonic models (e.g., Costamante et al. 2018), seem to challenge the scenario in which turbulence energizing the pre-accelerated electrons is the result of global instabilities excited by recollimation, since these instabilities are expected to strongly affect the flow, eventually disrupting it. For instance, in the simulations of Matsumoto et al. (2021) and Costa et al. (2024), the turbulent post-shock flow is strongly decelerated and the average bulk Lorentz factor is only a few, implying a moderate beaming of the radiation. Another factor that potentially impacts our scenario was revealed in recent results from the Imaging X-ray Polarimetry Explorer (IXPE). Observations of 1ES 0229+200 (Ehlert et al. 2023) show that the X-ray polarization is rather high, ∼18%, and strongly chromatic (i.e., X-ray polarization is higher than optical polarization). In our previous work, we suggested that in the model it is natural to expect a low polarization in all bands, but more accurate estimates are necessary in light of these detailed measurements. Following the phenomenological approach of Marscher & Jorstad (2022), the mean polarization degree of a turbulent emitting region is equal to

$$\begin{array}{c}\hfill ⟨P⟩\approx 0.75{[f_{\mathrm{ord}}^2+{(1-f_{\mathrm{ord}})}^2{N}_{\mathrm{cell}}^{-1}]}^{\frac{1}{2}},\end{array}$$(12)

where f_ord = B_ord/B (with B_ord the ordered component of the magnetic field) and N_cell is the number of turbulent cells. Fixing f_ord, the minimum is reached when N_cell → ∞; therefore, the maximum value of f_ord necessary to obtain the polarization degree observed by IXPE is given by f_ord, max = ⟨P⟩/0.75 ≈ 0.24, which implies δB²/B² ≳ 0.6, a value beyond the scope of the quasi-linear theory at the base of our treatment of stochastic acceleration.

In view of these difficulties, we propose a modification of the scenario that can potentially account for both the moderate relativistic Doppler factor and the chromatic polarization. While in the original idea the post-shock turbulence is supposed to be related to the onset of global jet instabilities triggered by recollimation, an alternative is to assume the existence of inhomogeneities in the upstream flow. This, in particular if the plasma is characterized by a low magnetization, promotes the development of Richtmyer–Meshkov-like instability at the shock, with the resulting onset of turbulence in the downstream region. Importantly, dedicated magnetohydrodynamic simulations show that close to the shock front the evolving turbulent eddies are small, but their size grows moving away from the shock (e.g., Mizuno et al. 2014). Since the jet does not mix with the external medium, it can maintain a moderate relativistic Doppler factor, as required by our fits. Moreover, the structure of the turbulence can explain chromatism: low-energy electrons that emit in the optical band are accelerated by small eddies (i.e., eddies characterized by a large k), which are present both near the shock front and far from it since large eddies naturally cascade at shorter lengths. On the other hand, high-energy electrons emitting in the X-ray band have large Larmor radii and are associated exclusively with the eddies with large wavelengths (small k) far from the shock front. We therefore expect that the effective volume where X-ray emission occurs is much smaller than that associated with the optical band, implying an X-ray polarization higher than the optical one, in line with observations (Marscher & Jorstad 2022). An improvement to our model could be the introduction of a time-dependent injection, to model the increase in the eddy size away from the shock front. Simulations are needed to understand the complete phenomenology of extreme TeV BL Lacs. Thus, our next step will be to use the Lagrangian particle framework to investigate these sources (Vaidya et al. 2018; Mukherjee et al. 2021).

The next generation of gamma-ray facilities, such as CTA and ASTRI, will provide more constraining measurements that will further reduce the model uncertainty. In particular, our model predicts a strong cutoff above 10 TeV that, if confirmed, will strongly limit the use of extreme TeV BL Lacs as beacons to test physics beyond the standard model (see, e.g., Galanti et al. 2020), such as axion-like particles and Lorentz invariance violation.

The exploration of the parameter space for phenomenological models has recently improved thanks to neural networks. After being trained on a sample of several SEDs, neural networks are capable of computing spectra in milliseconds, instead of the few seconds required by leptonic codes. Together with MCMC or nested sampling, neural networks can be used to efficiently explore the parameter space (Bégué et al. 2024; Tzavellas et al. 2024). We aim to implement this technique for our model as well.

Acknowledgments

We thank A. Mignone and S. Kundu for useful discussions on the numerical schemes. We acknowledge financial support by a INAF Theory Grant 2022 (PI F. Tavecchio) and the PRIN 2022 (2022C9TNNX) project. This work has been funded by the EU – Next Generation EU.

References

Appendix A: Corner plots

Fig. A.1. Corner plot of 1ES 0229+200.

Fig. A.2. Corner plot of 1ES 0347-121.

Fig. A.3. Corner plot of 1ES 1101-232.

Fig. A.4. Corner plot of RBG J0710+591.

Appendix B: MCMC of 1ES 1101-232 with δ = 40

Fig. B.1. Alternative corner plot of 1ES 1101-232.

Fig. B.2. Flux points with their errors (black) of 1ES 1101-232 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.

Appendix C: Comparison of numerical schemes

In previous papers, we used the Chang-Cooper (CC) scheme (Chang & Cooper 1970; Park & Petrosian 1996) to solve both equations of the system (1). We modified the boundary conditions for non-thermal electrons to no flux conditions. This update is negligible for the resulting emission, but it is useful to compare the performances of the CC algorithm with alternative schemes, such as the second-order implicit-explicit Runge-Kutta (IMEX-RK) algorithms (Kundu et al. 2021). Specifically, we tested the strong stability preserving scheme (2,2,2) of Pareschi & Russo (2005), but it presents several drawbacks for the scenario we are considering. IMEX-RK can be used exclusively for linear equations, so not for the turbulence equation, making the scheme only first-order accurate. Instead, CC can be adapted for nonlinear equations (Larsen et al. 1985). Therefore, we are stuck to the same grid accuracy, but, since the advection term is treated explicitly in IMEX-RK, several time steps are necessary, which makes the scheme much slower. However, we used IMEX-RK to test our CC implementation: we used the same momentum and wavenumber grid (i.e., 20 points per decade) but different time stepping (i.e. ,100 time steps for CC and C = 0.9 for IMEX-RK, where C is the Courant number). The two schemes produced comparable results (see Fig. C.1), confirming the proper functioning of CC implementation; therefore, we decided to keep employing it.

Fig. C.1. Comparison between two realizations of electron spectrum by IMEX-RK (solid red) and CC (solid blue) using the same parameters, momentum, and wavenumber stepping, but with shorter time stepping for the former.

All Tables

All Figures

	Fig. 1. Flux points with their errors (black) of 1ES 0229+200 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.
In the text

	Fig. 2. Flux points with their errors (black) of 1ES 0347−121 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.
In the text

	Fig. 3. Flux points with their errors (black) of 1ES 1101−232 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.
In the text

	Fig. 4. Flux points with their errors (black) of RGB J0710+591 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.
In the text

	Fig. 5. Final electron energy distribution for 1ES 0229+200 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior.
In the text

	Fig. 6. Final spectrum of the turbulence for 1ES 0229+200 with 1σ credible intervals (violet) and the median (blue) obtained from the model posterior.
In the text

	Fig. 7. Final electron timescales (acceleration, cooling, escape) for 1ES 0229+200 with 1σ credible intervals and the median obtained from the model posterior.
In the text

	Fig. 8. Final turbulence timescales (cascading and damping) for 1ES 0229+200 with 1σ credible intervals and the median obtained from the model posterior.
In the text

	Fig. 9. Distribution of the magnetization of 1ES 0229+200 obtained from the model posterior.
In the text

	Fig. A.1. Corner plot of 1ES 0229+200.
In the text

	Fig. A.2. Corner plot of 1ES 0347-121.
In the text

	Fig. A.3. Corner plot of 1ES 1101-232.
In the text

	Fig. A.4. Corner plot of RBG J0710+591.
In the text

	Fig. B.1. Alternative corner plot of 1ES 1101-232.
In the text

	Fig. B.2. Flux points with their errors (black) of 1ES 1101-232 with 1σ credible intervals (orange) and the median (red) obtained from the model posterior. Gray points correspond to Swift/UVOT data.
In the text

	Fig. C.1. Comparison between two realizations of electron spectrum by IMEX-RK (solid red) and CC (solid blue) using the same parameters, momentum, and wavenumber stepping, but with shorter time stepping for the former.
In the text