© A. Zech and M. Lemoine 2021

1. Introduction

Blazars – active galactic nuclei (AGN) with a jetted relativistic outflow pointing towards us – are commonly classified according to their characteristic double-humped spectral energy distribution (SED; Urry & Padovani 1995; Fossati et al. 1998). At the tip of the so-called blazar sequence (Padovani & Giommi 1995; Ghisellini et al. 2017), extreme blazars form a new population of BL Lacs (Costamante et al. 2001), designated as ‘ultra-high-frequency-peaked BL Lac objects’ (UHBLs) given the very high frequencies of their two emission peaks (Biteau et al. 2020). The sub-class of ‘extreme-TeV’ sources presents rather unusual properties, with the higher energy peak beyond ∼1 TeV and a hard spectrum in the γ-ray range. These objects still defy a theoretical interpretation in terms of the commonly used one-zone radiative models, which are generally successful in the representation of SEDs from BL Lac objects.

Although ten or so extreme-TeV sources are known to date (Biteau et al. 2020), current estimates based on data from the Fermi-LAT γ-ray telescope predict about five times more candidate sources that need to be confirmed with future instruments in the TeV range (Costamante 2020). For five of the known extreme-TeV blazars, largely simultaneous SEDs with a good multi-wavelength coverage from the infrared to the TeV band have become available, and they have already provided important clues on the physics of dissipation in these sources (Costamante et al. 2018).

Arriving at a satisfactory description of extreme-TeV sources in the standard one-zone leptonic framework seemingly requires two essential ingredients (e.g., Tavecchio et al. 1998; Ghisellini et al. 2002; Katarzyński et al. 2006; Costamante et al. 2018): (1) lower magnetic fields (≲10 mG) than for more common high-frequency-peaked BL Lac objects (HBLs), corresponding to magnetic field energy densities well below equipartition, and (2) a peculiar lepton distribution with most of the energy carried by particles of a large Lorentz factor, ${\overline{\gamma }}_{\mathrm{e}}\sim {10}^3-{10}^4$. The former requirement avoids an excessive softening of the γ-ray spectrum by synchrotron cooling and is necessary to reproduce the large separation between the synchrotron and synchrotron self-Compton (SSC) peak frequencies, while the latter is needed to accommodate the steep γ-ray spectrum and high peak frequency. Alternative models have been devised (e.g., Aharonian 2000; Böttcher et al. 2008; Lefa et al. 2011; Tavecchio 2014; Cerruti et al. 2015; Chhotray et al. 2017; Tavecchio & Sobacchi 2020), but they all come at the price of greater complexity and most often invoke other extreme parameter values, which often lack experimental or physical justification. For instance, lepto-hadronic models can provide reasonable SEDs with less extreme values of ${\overline{\gamma }}_{\mathrm{e}}$, but at the expense of a high jet power, mostly carried by the accompanying population of ultra-relativistic protons or the magnetic field (e.g., Cerruti et al. 2015).

In the present study we take at face value the low magnetisation inferred from the SEDs of extreme blazars and propose revisiting the one-zone model, adopting well-motivated microphysical prescriptions for particle acceleration, under the assumption that electrons are co-accelerated with protons in relativistic internal or recollimation shocks. Two key motivations of our study are the following: (1) at low magnetisation, relativistic shock acceleration is expected to dominate the physics of dissipation and of the generation of non-thermal power laws, and (2) electrons can be efficiently preheated in the transition layer of electron-ion relativistic shocks, up to a fraction of the thermal energy of shocked ions, implying large minimum Lorentz factors for the electron population, as seemingly required by phenomenology (e.g., Vanthieghem et al. 2020, and references therein). We further aim at making our approach as self-contained as possible, by fixing the microphysical parameters describing the acceleration process from first principle studies of relativistic shock models. Our main objective is to provide a critical discussion of such models to reproduce the SEDs of extreme blazars, in a modern context. This contrasts with generic one-zone models, which instead parameterise the accelerated particle population and adjust the parameters to the observed SEDs.

We organise the discussion as follows. In Sect. 2 we discuss existing constraints on the microphysical parameters of shock acceleration. The standard one-zone SSC model is then extended to include a population of protons, and physical constraints from the considered acceleration mechanism are translated into constraints on the parameters that describe the two particle populations. Section 3 then applies the resulting model to the set of well-covered SEDs. In a first scenario, generic solutions with a minimum number of free parameters are presented. A second scenario applies realistic parameters for acceleration on a mildly relativistic internal or recollimation shock, while a third scenario accounts for the possibility of re-acceleration at multiple shock waves, which appears necessary to reproduce the hard γ-ray spectra of certain objects. A critical discussion of this model, as well as the direct consequences for our understanding of the jets of extreme blazars, is given in Sect. 4.

2. A theoretical model for e-p co-acceleration in blazar internal shocks

2.1. A large ${\overline{\gamma }}_{\mathrm{e}}$

We define the Lorentz factor ${\overline{\gamma }}_{\mathrm{e}}$ as the average over the population,

$$\begin{array}{c}\hfill {\overline{\gamma }}_{\mathrm{e}}=\frac{1}{n_{\mathrm{e}}}{\displaystyle \int }\mathrm{d}{\gamma }_{\mathrm{e}}{\gamma }_{\mathrm{e}}\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}{\gamma }_{\mathrm{e}}},\end{array}$$(1)

with n_e = ∫dγ_e dn_e/dγ_e the electron number density. Blazar SEDs are generically well represented by radiative leptonic emission from a (broken) power-law population with spectral indices s₁ (resp. s₂) below (resp. above) a break Lorentz factor γ_e, b. For such a distribution of Lorentz factors, ${\overline{\gamma }}_{\mathrm{e}}\sim {\gamma }_{\mathrm{e},\min }{({\gamma }_{\mathrm{e},\mathrm{b}}/{\gamma }_{\mathrm{e},\min })}^{2-s_1}$ if s₁ < 2 < s₂, and ${\overline{\gamma }}_{\mathrm{e}}\sim {\gamma }_{\mathrm{e},\min }$ if 2 < s₁ < s₂. The energy density u_e of the population can then be written $u_{\mathrm{e}}\sim {\overline{\gamma }}_{\mathrm{e}}{n}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2$.

Therefore, to convert a population of electrons (or pairs) of initial Lorentz factor γ₀ to the above power-law spectrum necessitates the dissipation of an energy reservoir that is ${\overline{\gamma }}_{\mathrm{e}}/{\gamma }_0$ times larger than the initial particle energy. In this sense, the large inferred value of ${\overline{\gamma }}_{\mathrm{e}}$ provides an important clue regarding the amount of energy that is dissipated and its origin, given that, in the models of Costamante et al. (2018), ${\overline{\gamma }}_{\mathrm{e}}\sim {10}^3-{10}^4$ arises as a generic prediction independently of the value of s₁.

In a (co-moving) magnetic field of strength B = 10 B₋₂ mG, electrons of Lorentz factor γ_e = 10⁴ γ_e 4 cool through synchrotron radiation on a length scale

$$\begin{array}{c}\hfill l_{\mathrm{c}}\simeq 8{\mathrm{\Gamma }}_{\mathrm{j}}{B}_{-2}^{-2}{\gamma }_{\mathrm{e}4}^{-1}\mathrm{pc},\end{array}$$(2)

as measured in the source rest frame, and expressed in terms of Γ_j the jet Lorentz factor. Hence, one cannot a priori exclude that electrons gained their large average Lorentz factors in some primary dissipation region, located well upstream of the acceleration region where the power-law population is generated. Such a scenario is nevertheless constrained by the apparent lack of radiation associated with that putative primary dissipation process, as well as by the shortened cooling length in the primary dissipation region with enhanced magnetic strength. In the absence of a definite model for the evolution of the magnetic field strength, electron distribution etc. along the jet, it is difficult to make this statement more quantitative as the synchrotron luminosity scales with the combination ${\mathrm{\Gamma }}_{\mathrm{j}}^4{n}_{\mathrm{e}}{R}^3{u}_B{\gamma }_{\mathrm{e},\min }^{s-1}{\gamma }_{\mathrm{e},\max }^{3-s}$ for a single power law of index s < 3 extending from γ_e, min to γ_e, max, where u_B denotes the magnetic energy density, n_e the electron density, R the jet radius, and Γ_j its Lorentz factor. We thus leave this possibility open in the following.

In this context, it is particularly interesting to investigate the alternative possibility that the heating of electrons to large minimal Lorentz factors is, in fact, associated with the acceleration process itself. In the following, we focus on the case of shock acceleration, while alternative scenarios involving reconnection or turbulent acceleration are examined in Appendix A.

2.2. Electron heating in relativistic shocks

Consider a shock front, moving at Lorentz factor γ_sh (and velocity β_shc) with respect to the un-shocked plasma. In the context of extreme blazars, we are mostly interested in the case of mildly relativistic or truly relativistic shocks (i.e. u_sh ≡ γ_shβ_sh ≳ 1). We also assume here a shock of normal incidence, meaning that the upstream flow enters the shock front along the shock normal, whose direction is perpendicular to the shock front. We generalise this to the case of oblique shocks in Sect. 2.4.2.

If the composition is purely leptonic, the electrons are energised through shock crossing according to the standard Rankine-Hugoniot conditions, or rather their relativistic generalisation (e.g., Kirk & Duffy 1999, and references therein). For a strong, sub-relativistic pair shock, the post-shock temperature reads $k_{\text{B}}{T}_{\text{e}}/m_{\text{e}}{c}^2\simeq 0.19u_{\text{sh}}^2$, while for a strong, relativistic pair shock, k_BT_e/m_ec² ≃ 0.24u_sh. For reference, k_BT_e/m_ec² ≃ 0.7 for u_sh = 3 in the mildly relativistic regime, close to the fully relativistic limit.

This situation changes rather dramatically in an electron-ion shock, assuming for simplicity an equal number of ions and electrons. Most of the energy that enters the shock front is now carried by the ions, but a fraction of it is given to the electrons, which are thereby preheated up to a fraction of equipartition. More specifically, in the reference frame in which the shock front lies at rest, the kinetic energy density in (cold) protons is e_p = γ_sh(γ_sh−1)n_pm_pc² (n_p is understood as a proper density here), corresponding to an energy-per-particle ⟨E_p⟩=(γ_sh−1)m_pc², and similarly for the electrons with _p ↔ _e. Efficient preheating of the electrons means that ⟨E_e⟩ becomes a substantial fraction of ⟨E_p⟩, instead of ⟨E_e⟩/⟨E_p⟩=m_e/m_p ≪ 1, which would be expected if both species were to satisfy the fluid shock crossing conditions independently from each other.

The mechanism(s) through which electrons gain energy at the expense of ions in shock fronts remain somewhat elusive. In the relativistic, weakly magnetised limit, one promising scenario is that in which electrons undergo collisionless Joule heating, driven by a longitudinal electric field and by the effective gravity that results from the slowdown of the microturbulence that is self-generated in the shock precursor (see e.g., Lemoine et al. 2019c; Vanthieghem et al., in prep.). Phenomenological models of gamma-ray burst afterglows provide a nice empirical confirmation of this energy transfer between ions and electrons, as the electron energy fraction is almost always found to be within a factor of a few from the ion energy fraction (e.g., Kumar & Zhang 2015, and references therein).

In the absence of a definitive model, we can rely on particle-in-cell (PIC) simulations to provide estimates for the average electron Lorentz factor. Generally speaking, the preheating appears more efficient at lower magnetisations σ (Sironi et al. 2013), likely because of a larger precursor and because of the presence of intense Weibel-like microturbulence. For reference, we define σ as the ratio of the magnetic energy density to the particle energy density in the rest frame of the un-shocked (ambient) plasma. Simulations also report a more efficient preheating in the fully relativistic regime – meaning γ_sh ≳ 3 − 5 – than in the mildly relativistic regime γ_sh ∼ 1 − 3. For example, Crumley et al. (2019) measure T_e ≃ 0.3T_p and k_BT_p ≃ 0.2m_pc² for a mildly relativistic sub-luminal shock of velocity β_sh = 0.8, corresponding to γ_sh = 1.7, at magnetisation σ ∼ 0.007 (here defined downstream). This electronic temperature corresponds to a mean Lorentz factor ${\overline{\gamma }}_{\mathrm{e}}\sim 300$. At larger shock velocities, and smaller magnetisations, Sironi & Spitkovsky (2011) and Sironi et al. (2013) find that the fraction of energy stored in electrons increases up to half of that in the ions, implying k_BT_e ∼ 0.1γ_shm_pc², or a mean Lorentz factor ${\overline{\gamma }}_{\mathrm{e}}\sim 600{\gamma }_{\mathrm{sh}}$.

In the sub-relativistic limit, and more generally in cases in which electron preheating is not efficient, the electrons suffer from a well-known injection problem into the Fermi acceleration cycle. In such cases, the energy density of the suprathermal tail of accelerated electrons represents a modest fraction ϵ_e of the energy that is incoming into the shock, for example ϵ_e ≃ 0.001 for the case with β_sh = 0.8 discussed in Crumley et al. (2019). Correspondingly, the minimum Lorentz factor of the suprathermal electron power law – that is, the Lorentz factor at which the power law emerges out of the thermal Maxwellian – is substantially larger than ${\overline{\gamma }}_{\mathrm{e}}$. In the relativistic limit, ϵ_e ≳ 0.1 because preheating, hence injection is efficient and the minimum Lorentz factor of the accelerated power law is a factor of a few larger than ${\overline{\gamma }}_{\mathrm{e}}$, at most.

If electron preheating were inefficient, the electron distribution would be dominated by the thermal bump, hence it would not account satisfactorily for the observed SEDs of extreme blazars, which not only require a large ${\overline{\gamma }}_{\mathrm{e}}$, but also a hard spectrum with index s ∼ 2. In the following, we thus adopt ϵ_e = 0.1 and ${\overline{\gamma }}_{\mathrm{e}}\sim 600{\gamma }_{\mathrm{sh}}$ as fiducial values, corresponding to relativistic, weakly magnetised electron-ion shocks. We also assume one electron per ion, as larger electron multiplicities are expected to lead to lower mean energy per particle, although this case has not been explicitly probed by numerical experiments.

2.3. Fixing parameters from the microphysics of shock acceleration

In the fully relativistic regime, meaning a shock Lorentz factor γ_sh ≳ 3 − 5, particle acceleration appears restricted to the weakly magnetised limit, σ ≲ 10⁻³ (see e.g., Vanthieghem et al. 2020 and references therein). Strictly speaking, this result applies to super-luminal shocks, but highly relativistic shocks are generically super-luminal as a consequence of the Lorentz boost (a factor γ_sh) of the perpendicular magnetic field components in the shock rest frame. In the mildly relativistic limit, however, such effects disappear and sub-luminal configurations become about as likely. Consequently, particle acceleration can become efficient, even at relatively high magnetisations (see e.g., Crumley et al. 2019).

In the regime of low magnetisation that we are interested in, the size of the precursor is governed by the scattering of particles in the electromagnetic micro-turbulence that they themselves excite through beam-plasma instabilities (e.g., Plotnikov et al. 2018), and which is responsible for the shock transition itself. This general problem has been studied in some detail in the limit of negligible initial magnetisation (e.g., Lemoine et al. 2019a,b,c; Pelletier et al. 2019). The microturbulence is mostly excited by the current filamentation instability (often termed Weibel), which generates an intense magnetic field on skin depth scales λ_p, with ${\lambda }_{\mathrm{p}}\equiv c{(4\pi n_{\mathrm{p}}{e}^2/m_{\mathrm{p}})}^{-1/2}\sim {10}^7{n}_{p0}^{1/2}$cm (n_p 0 = n_p/1 cm⁻³). The effective magnetisation carried by this turbulence, written ϵ_B, generally saturates at values ϵ_B ∼ 0.01 within thousands of skin depths from the shock front. Far downstream from the shock, it is expected to decay with distance as a mild power law (Lemoine 2013) through collisionless damping. Consequently, the effective magnetisation in the radiation region, which we write σ_rad, is expected to be smaller than 10⁻², but at least as large as the initial σ value.

The modelling of gamma-ray burst afterglows offers an interesting perspective on the magnitude of σ_rad. In those sources, the external magnetisation is significantly lower (σ ∼ 10⁻⁹) than expected here, yet σ_rad is inferred to be as large as ∼10⁻⁵ − 10⁻⁴ (Lemoine et al. 2013; Santana et al. 2014). That effective magnetisation is generally understood as a left-over of the shock-generated field, or as a the result of further buildup through additional instabilities, for example the Richtmyer–Meshkov instability at the shock front (Inoue et al. 2011) or a Rayleigh-Taylor instability at the contact discontinuity (Duffell & MacFadyen 2014). The value of σ_rad that we use in our models is comparable to that above; in that sense, gamma-ray burst modelling provides empirical support to our model.

At moderate magnetisation, 10⁻⁵ ≲ σ ≲ 10⁻², a different current-driven instability is expected to contribute to the generation of the magnetised microturbulence (Lemoine et al. 2014). The consequences of this instability, which is driven by the perpendicular current carried by the accelerated particles around the mean background magnetic field, have not been examined to the extent needed. In a first approximation, we can consider that the turbulence will have the same general characteristics as that generated by the current filamentation instability, that is, a strength ϵ_B ∼ 0.01 and a length scale ∼ 𝒪(λ_p), and we do so in the following.

The properties of the microturbulence in the vicinity of the shock govern the acceleration process, hence the parameters that characterise the spectrum of accelerated particles. In particular, the spectral index is expected to take values s ≃ 2 − 2.3. For definiteness, we assume s ≃ 2.2 in our application to extreme-TeV blazars, corresponding to mildly relativistic shocks. The minimum Lorentz factor is ${\gamma }_{\mathrm{e},\min }\simeq {\overline{\gamma }}_{\mathrm{e}}\simeq 600{\gamma }_{\mathrm{sh}}$ as discussed before and the maximal Lorentz factor γ_e, max is determined by the conjunction of a number of constraints: (1) the age limit, t_acc ≲ d/(Γ_jβ_jc), where t_acc represents the acceleration timescale to γ_e, max and d the distance from the jet base, Γ_j (resp. β_j) the jet Lorentz factor (resp. velocity); (2) the radiative loss limit, t_acc ≲ t_syn, written here as a competition with the synchrotron loss timescale; (3) the scattering limit (for super-luminal shocks) t_scatt ≲ r_L, 0, where r_L, 0 represents the particle gyroradius in the coherent (background) magnetic field. We note that the limit associated with lateral escape through the jet boundary, viz.t_acc ≲ r, is comparable to the age limit, if r ∼ Θ_jd and the jet opening angle Θ_j ∼ 1/Γ_j.

The acceleration timescale is given by t_acc ≃ t_scatt in relativistic shocks, where t_scatt denotes the scattering timescale, which departs from the Bohm regime (t_scatt ∝ r_L) at weak magnetisation (Vanthieghem et al. 2020, and references therein). Consequently, the radiative limit can be written γ_e,max ∼ 10⁻⁷(n_e/1 cm⁻³)^−1/6, which does not depend on the strength of the magnetic field in the acceleration region, contrary to the usual relation γ_e, max ∝ B^−1/2 obtained for Bohm scaling.

The scattering constraint ensures that particles are able to scatter effectively in the microturbulence, hence to cross the shock repeatedly, before being advected away from the shock by their gyration around the regular magnetic field lines (Pelletier et al. 2009; Lemoine & Pelletier 2010). This constraint thus assumes the magnetic field to be super-luminal and it would not arise in sub-luminal configurations. In this sense, it can be regarded as conservative. For scattering in microturbulence, this constraint can be reexpressed as a function of magnetisation:

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{e},\max }\lesssim \frac{{\gamma }_{\mathrm{e},\min }}{\sqrt{\sigma }}.\end{array}$$(3)

Its main effect is to reduce the dynamic range over which particles can be accelerated. Its dependence on σ derives from the scalings t_scatt ∝ (γ_e/γ_sh)² and r_L, 0 ∝ (γ_e/γ_sh)/⟨B⟩. The maximum Lorentz factor arising from this limit is multiplied by m_p/m_e in an electron-proton plasma. (See Vanthieghem et al. 2020 for more details.)

In the models that we adjust to the blazar SEDs further below, we find that this third constraint provides the limiting factor, although the radiative constraint and the age constraints do not lag by much. To put in some numbers, assume ϵ_B = 10⁻², σ = 10⁻⁵, γ_sh = 3, n_e = 1 cm⁻³. Then, the scattering constraint gives γ_e, max ∼ 2 × 10⁵, while the radiative and age constraints both lead to γ_e, max ∼ 10⁷. The severity of the scattering constraint results from the relatively “low” value of γ_e, min, here 1800; below, we discuss models in which γ_e, min is larger because of acceleration at multiple shocks. In such cases, the maximal Lorentz factor that results from that scattering constraint becomes large enough to explain radiation up to TeV energies. Again, we stress that this scattering limit is conservative as it formally applies to super-luminal shocks and not sub-luminal ones.

We also modelled the radiation from protons, although we find that it is always negligible for the model parameters that we derive. The spectrum characterising the protons can be similarly defined by a minimum Lorentz factor, γ_p, min ∼ γ_sh, a maximal Lorentz factor determined by the conjunction of the age and scattering constraints above and the same spectral index as for the electrons.

2.4. Shock origin and geometry

2.4.1. Shocks of normal incidence

The interaction of a jet with an obstacle generally triggers a double shock configuration, including a forward and a reverse shock. Depending on the relative momentum fluxes carried by the jet and the obstacle, either or both of these shocks can be strong or weak (Sari & Piran 1995). Given that particles release their radiation in the downstream rest frame of the shock where they are accelerated, the velocity of that downstream frame sets the Doppler factor that modulates the radiation. For extreme-TeV blazars, the large inferred Doppler boosting indicates that this downstream frame moves at large relativistic velocities towards the observer in the source frame. For shocks of normal incidence, this implies in turn that the shock itself moves at large velocities in the source frame, because the downstream rest frame moves at sub-relativistic velocities relative the shock.

Consequently, if the jet overtakes an obstacle, we need to assume that the jet ram pressure far exceeds that of the obstacle, so that the Doppler boosting of the emission is given by that of the jet itself. The shock Lorentz factor γ_sh corresponds here to the relative Lorentz factor between the shock front and the obstacle, and it is of the order of Γ_j if the obstacle moves at slow velocities in the source rest frame. This may give rise to shocks of large Lorentz factors. However, obstacles can penetrate a jet only if their ram pressure exceeds that of the jet (e.g., Araudo et al. 2010; Barkov et al. 2010; Bosch-Ramon et al. 2012), which would instead imply a weak forward shock.

Alternatively, if an ‘obstacle’ overtakes the jet, as in the ‘blob-in-jet’ scenario, the Doppler factor will be at least as large as that of the jet itself. If the blob ram pressure exceeds that of the jet, the Lorentz factor of the emission region is given by that of the blob, Γ_b, while the shock Lorentz factor γ_sh is set as before by the relative Lorentz factor between the blob and the jet. From γ_sh = Γ_jΓ_b(1 − β_jβ_b), Γ_j ≫ 1, Γ_b ≫ 1 and Γ_b ≫ Γ_j, it follows

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{sh}}\sim {\mathrm{\Gamma }}_{\mathrm{b}}/(2{\mathrm{\Gamma }}_{\mathrm{j}}).\end{array}$$(4)

Mildly relativistic shocks can be attained for sufficiently large ratios of Γ_b to Γ_j. We consider such a configuration in Sect. 3.2.

2.4.2. Oblique shocks

In the source rest frame, the shock(s) may well be of oblique incidence, especially if they are recollimation shocks (e.g., Perucho 2013; Martí et al. 2016, and references therein). The general properties of these recollimation shocks have been explored in various studies (e.g., Dubal & Pantano 1993; Komissarov & Falle 1997; Nalewajko & Sikora 2009; Bromberg & Levinson 2009). We summarise in Appendix B the characteristics of interest to us, in particular the jump conditions, and we establish here the correspondence between these properties and the above constraints on particle acceleration.

In brief, these constraints remain applicable to oblique shocks once they are expressed in the so-called shock normal rest frame, in which the (un-shocked) flow impinges on the shock at normal incidence (Begelman & Kirk 1990). The Lorentz factor γ_sh, which corresponds to the Lorentz factor of the shock with respect to the upstream plasma in this shock normal frame, is written Γ_j < |n in Appendix B. It can be related to the relative Lorentz factor of the (pre-shock) jet with the (oblique) shock surface, written Γ_j<, through Eq. (B.2): γ_sh = Γ_j</Γ_n|s, where Γ_n|s represents the Lorentz factor associated with the boost to the shock normal frame. The main result here can be phrased in terms of the angle |θ_j<−α| between the flow incidence and the shock surface (θ_j< and α, respectively, represent the angles of the flow and the shock surface with respect to the jet axis). Given that the jet opening angle is of the order of 1/Γ_j, one expects |θ_j<−α| ∼ 1/Γ_j in a simple re-confinement scenario (Nalewajko & Sikora 2009), and we therefore write |θ_j<−α| ∼ κ/Γ_j<. Then, we obtain ${\gamma }_{\mathrm{sh}}\simeq \sqrt{1+{\kappa }^2}$. The shock is thus mildly relativistic for κ of the order of unity. The angle |θ_j<−α| may however show a different behaviour in jets where recollimation shocks arise from a difference in pressure between jet components and the ambient medium (e.g., Fichet de Clairfontaine et al. 2021).

Similarly, we can write the Lorentz factor of the flow past the shock as ${\mathrm{\Gamma }}_{\mathrm{j}>}\simeq {\mathrm{\Gamma }}_{\mathrm{n}|\mathrm{s}}\simeq {\mathrm{\Gamma }}_{\mathrm{j}<}/\sqrt{1+{\kappa }^2}$. This Lorentz factor is that which controls the amount of Doppler boosting from the emission region, provided the recollimation shock is stationary in the source frame. We note that the flow is refracted by an angle of the order of |θ_j<−α| in the source frame as it crosses the shock, which may therefore affect the overall Doppler boosting from the post-shock region quite substantially relative to that in the pre-shock region. We do not consider this effect in the models that follow, for the sake of simplicity.

In our models below, we seek to describe the radiation emitted by an element of plasma crossing such a shock, or a series of shocks. In this description, we can model the plasma element as a blob containing a fixed number of particles, in which particles get accelerated at some stage then radiate. Even though we are dealing with a standing shock, which provides a steady emission pattern, the corresponding Doppler boosting of the radiation from the blob, as seen by the distant observer, is a factor δ⁴ (cf. Sikora et al. 1997), with δ the bulk Doppler factor. For the sake of completeness, we provide more details on this issue in Appendix C.

2.4.3. Interaction with multiple shocks

As recollimation (or, more generally, standing) shocks may come in series of shock fronts in the jets of radio-galaxies, we must consider the possibility that the radiating blob undergoes multiple episodes of shock acceleration, which may substantially modify the spectrum of accelerated particles (White 1985; Achterberg 1990; Schneider 1993; Pope & Melrose 1994; Meli & Biermann 2013). Particles are heated and accelerated in those shocks, but they also suffer adiabatic cooling (and possibly, radiative cooling) in the rarefaction regions that separate them.

As a particle crosses a shock front, its energy in the local frame of rest of the plasma is increased on average by the relative Lorentz factor between pre-shock and post-shock flows. Average is understood here over the orientation of the particle momentum in the pre-shock flow. For an oblique shock, the relative Lorentz factor Γ_rel can be written in terms of normal shock frame quantities Γ_rel ≃ 0.7Γ_j < |n, where the prefactor holds for mildly relativistic shocks with Γ_j < |n ≳ 2 − 3. In terms of source frame quantities, Γ_rel ≃ 0.5Γ_j</Γ_j> (in analogy with Eq. (4)).

Adiabatic cooling in the rarefaction regions implies, for one-dimensional dilatation along the flow axis, $\mathrm{\Delta }lnp=-\frac{1}{3}\mathrm{\Delta }ln{u}_{\mathrm{j}}$, with u_j the flow four-velocity. Consider for instance a velocity profile such that Γ_j< decreases abruptly to Γ_j> through a recollimation shock, then steadily re-increases to its initial value Γ_j< through rarefaction. Overall, the momentum of the particle increases by $\mathit{g}\sim {\mathrm{\Gamma }}_{\mathrm{rel}}^{2/3}$ between a time defined as immediately before the crossing of the first shock and that immediately before the crossing of the second shock. The difference with the sub-relativistic regime studied in the above references should be noted : there the particle gains little energy through shock crossing, but cools efficiently in the rarefaction region, leading to an overall energy loss in between two consecutive shock crossings.

Beyond this sequence of heating and cooling, which shifts the energy spectrum in momentum space while preserving its shape, the particles are also re-accelerated at the shock. This now modifies the spectral shape. This effect can be described, in a first approximation, through the effective propagator

$$\begin{array}{c}\hfill G({\gamma }_{>},{\gamma }_{<})=\frac{s-1}{g{\gamma }_{<}}{\left(\frac{{\gamma }_{>}}{g{\gamma }_{<}}\right)}^{-s}\mathrm{\Theta }({\gamma }_{>}-g{\gamma }_{<}),\end{array}$$(5)

which provides the probability distribution of the outgoing particle Lorentz factor γ_> consequent to crossing a second shock, assuming a Lorentz factor γ_< immediately after crossing a first shock. Here g represents the energy gain from one shock to the next, $\mathit{g}\simeq {\mathrm{\Gamma }}_{\mathrm{rel}}^{2/3}$ as discussed above. The post-shock particle distribution can be obtained through the convolution of the pre-shock distribution with this propagator, which gives, after n shock crossings (White 1985; Achterberg 1990; Schneider 1993)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{>}^{(n)}}{\mathrm{d}{\gamma }_{>}}=\frac{{(s-1)}^{n+1}}{n!g^n{\gamma }_{\min }}{\left(\frac{{\gamma }_{>}}{g^n{\gamma }_{\min }}\right)}^{-s}ln{\left(\frac{{\gamma }_{>}}{g^n{\gamma }_{\min }}\right)}^n.\end{array}$$(6)

This has two important consequences: (1) the spectrum hardens because of reacceleration; (2) the effective injection Lorentz factor, which was γ_min at the first shock, has become gⁿγ_min at the nth shock. The hardening is stronger at momenta close to the effective injection momentum than at high energy (HE). In particular, the maximum of $\text{d}{N}_{>}^{(n)}/\text{d}\ln \gamma$, which provides a reasonable estimate for ${\overline{\gamma }}_{\mathrm{e}}$ is found at γ  =  gⁿ γ_minexp[n/(s−1)], which indicates a rapid rise of the effective ${\overline{\gamma }}_{\mathrm{e}}$ with the number of shock crossings, all the more so if g is larger than unity (i.e. for relativistic shocks).

To make connection with the model parameters, γ_e, min is now given by gⁿ γ_e, min, with γ_e, min ≃ 600γ_sh and the spectrum is given by Eq. (6) above.

2.5. One-zone model with e − p co-acceleration in relativistic shocks

The standard one-zone SSC model has usually nine free parameters: the strength of the uniform magnetic field B, the size of the spherical emission region R, its Doppler factor δ, the minimum, maximum and break Lorentz factors of the electron distribution γ_e, min, γ_e, max, γ_e, br, the two indices of the electron distribution before and after the break s₁ and s₂, and its normalisation K_e at γ = 1¹. Even with a well-sampled multi-wavelength dataset, this model remains highly degenerate.

Tavecchio et al. (1998) have studied the constraints on these parameters for a set of six observables extracted from the SED: the peak frequencies and luminosities of the synchrotron and inverse Compton (IC) components, ν_syn, ν_IC, ν_synL_syn(ν_syn), ν_ICL_IC(ν_IC), and the indices in the differential photon flux F_ν before and after the synchrotron peak, α₁ and α₂. A detailed method for an exploration of this parameter space for a given dataset, allowing for statistical uncertainties in the data, has been discussed by Cerruti et al. (2013).

The constraints on the magnetisation, on γ_e, min, γ_e, max and on the spectral index, which we derived above from the microphysics of relativistic shock acceleration, remove most of the inherent degeneracy of the standard model. We demonstrate this by considering a power-law electron distribution with an exponential cutoff of the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{e}}}{\mathrm{d}\gamma }=K_{\mathrm{e}}{\gamma }^{-s}{e}^{-\gamma /{\gamma }_{\mathrm{e},\max }}.\end{array}$$(7)

This simple parameterisation turns out to be well suited to characterise the SEDs of extreme-TeV blazars.

For an initially chosen value of δ, first the maximum electron Lorentz factor γ_e, max can be determined from the IC peak frequency ν_IC. In the Klein-Nishina regime, which is relevant for the datasets we want to investigate in the following section, this relation is given following Tavecchio et al. (1998):

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{e},\max }\simeq \frac{1}{g({\alpha }_1,{\alpha }_2)}\frac{1+z}{\delta }\frac{h{\nu }_{\mathrm{IC}}}{m_{\mathrm{e}}{c}^2},\end{array}$$(8)

with

$$\begin{array}{c}\hfill g({\alpha }_1,{\alpha }_2)=exp(\frac{1}{{\alpha }_1-1}+\frac{1}{2({\alpha }_2-{\alpha }_1)}).\end{array}$$(9)

Setting s ≃ 2.2 then implies α₁ ≃ 0.6 as α = (s − 1)/2 because the particles do not effectively cool through synchrotron radiation. The value of α₂ is fixed to 1.75 for all sources under study, which provides a sufficiently good approximation for the shape of the synchrotron spectrum just after the peak, when assuming an electron distribution with exponential cutoff.

With γ_e, max known, the magnetic field strength B in the radiation region can be estimated from the value of ν_syn, giving

$$\begin{array}{c}\hfill B=27\mathrm{mG}\frac{1+z}{\delta }\frac{{\nu }_{\mathrm{syn}}}{{10}^{17}\mathrm{Hz}}{\left(\frac{{\gamma }_{\mathrm{e},\max }}{{10}^6}\right)}^{-2}.\end{array}$$(10)

Since the above formulas were derived for a broken-power-law distribution, a small adjustment is needed for a good representation with our electron distribution. Multiplying the observed values of ν_syn and ν_IC by a factor 2.0 has been found to adjust well for the gradual flux decrease due to the exponential cutoff.

Provided the maximum electron energy is set by the scattering constraint, following Eq. (3), and the minimum electron energy can be estimated, at least broadly, from the low-frequency part of the SED and the spectral slope in the GeV range, the magnetisation of the plasma flow can be estimated as σ ≃ (γ_e, min/γ_e, max)².

The size of the homogeneous emission region is derived from the peak luminosities, following again Tavecchio et al. (1998) in the Klein-Nishina regime:

$$\begin{array}{c}\hfill R={\left\{2\frac{f({\alpha }_1,{\alpha }_2)}{{\delta }^4}\frac{{\left[{\nu }_{\mathrm{syn}}{L}_{{\nu }_{\mathrm{syn}}}\right]}^2}{{\nu }_{\mathrm{IC}}{L}_{{\nu }_{\mathrm{IC}}}{B}^{2}c}{\left(\frac{3\delta m_{\mathrm{e}}{c}^2}{4{\gamma }_{e,max}h{\nu }_{\mathrm{syn}}}\right)}^{(1-{\alpha }_1)}\right\}}^{1/2},\end{array}$$(11)

with

$$\begin{array}{c}\hfill f({\alpha }_1,{\alpha }_2)=\frac{1}{1-{\alpha }_1}+\frac{1}{{\alpha }_2-1}.\end{array}$$(12)

We may then verify our assumption of a single power-law distribution without cooling break, by comparing the frequency at which synchrotron cooling becomes relevant, namely

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{c}}\simeq 7.7\times {10}^8{\left(\frac{B}{1\mathrm{mG}}\right)}^{-2}{\left(\frac{R}{3\times {10}^{16}\mathrm{cm}}\right)}^{-1}\end{array}$$(13)

versus γ_e, max. Given the low magnetisation required for efficient shock acceleration, the assumption holds in our application. IC cooling is not significant either in the absence of strong Compton dominance.

Hence, the only missing parameter is the normalisation of the electron distribution K_e. Following our scenario outlined in the previous section, a population of protons is added to the model. In a standard one-zone lepto-hadronic model, this would lead to four additional free parameters for a power-law distribution (e.g., Cerruti et al. 2013; Zech et al. 2017). In our scenario, these parameters are constrained by the microphysics, as discussed earlier.

The normalisation K_e can then be determined from the magnetisation. We first determine the normalisation of the proton spectrum K_p from

$$\begin{array}{c}\hfill u_{\mathrm{p}}\equiv {\displaystyle {\int }_{{\gamma }_{\mathrm{p},\min }}^{{\gamma }_{\mathrm{p},\max }}\mathrm{d}\gamma K_{\mathrm{p}}{\gamma }^{-s+1}{e}^{-\gamma /{\gamma }_{\mathrm{p},\max }}{m}_{\mathrm{p}}{c}^2=\frac{u_B}{{\sigma }_{\mathrm{rad}}},}\end{array}$$(14)

with σ_rad the value of the magnetisation in the radiation region, which we recall is expected to be equal or larger than the pre-shock magnetisation. For a pure electron-proton plasma with equal particle number densities n_e = n_p, we then find

$$\begin{array}{c}\hfill K_{\mathrm{e}}=K_{\mathrm{p}}\frac{{\int }_{{\gamma }_{\mathrm{p},\min }}^{{\gamma }_{\mathrm{p},\max }}\mathrm{d}\gamma {\gamma }^{-s}{e}^{-\gamma /{\gamma }_{\mathrm{p},\max }}}{{\int }_{{\gamma }_{\mathrm{e},\min }}^{{\gamma }_{\mathrm{e},\max }}\mathrm{d}\gamma {\gamma }^{-s}{e}^{-\gamma /{\gamma }_{\mathrm{e},\max }}}.\end{array}$$(15)

The normalisation K_e does not a priori fit the absolute level of the observed flux defined by ν_synL_νsyn and ν_ICL_νIC. Thus, in an iterative procedure, the value of δ needs to be varied to search for a solution for a given SED. The minimum allowed value of δ can be determined from constraints on the source opacity and from observations on the variability timescale, while its maximum value is more loosely limited by energetics and by the statistics of beamed and un-beamed sources (Tavecchio et al. 1998).

In Eq. (14) we have assumed that the proton and electron spectra share a same spectral index. This is a rather common assumption, of no consequence here, because the hadronic contribution to the radiation will be found to be undetectable in the present model. For reference, we assume γ_p, min = γ_sh and ${\gamma }_{\mathrm{p},\max }={\gamma }_{\mathrm{p},\min }/\sqrt{\sigma }$, in line with the scalings adopted for the electrons, up to preheating in the shock transition. When considering a scenario where particles are re-accelerated on consecutive shocks, an analogous procedure can be followed by replacing the power-law particle distributions with the distribution given by Eq. (6).

To enable a comparison of the energy requirements of the co-acceleration model to that of standard leptonic models, the jet power is evaluated using the usual approximation (for a two-sided jet):

$$\begin{array}{c}\hfill L_{\mathrm{j}}\approx 2\pi R^2\beta c{\mathrm{\Gamma }}^2(u_B+u_{\mathrm{e}}+u_{\mathrm{p}}+u_{\gamma }),\end{array}$$(16)

where u_γ represents the energy density contained in the SSC radiation field.

We assume a value for the bulk Lorentz factor Γ = 0.5 δ, corresponding to a situation where the jet is closely aligned with the line of sight. The required luminosity L_j needs to be compared to the theoretically available power from the central engine, for which the Eddington luminosity of the central black hole can be used as an order-of-magnitude estimate.

3. Application to extreme-TeV blazars

The co-acceleration scenario was applied to five of the six multi-wavelength SEDs from extreme-TeV blazars presented by Costamante et al. (2018) for the sources 1ES 0229+200, 1ES 1101−232, 1ES 0347−121, RGB J0710+591, 1ES 1218+304, and 1ES 0414+009. These datasets were chosen due to the dense and homogeneous wavelength coverage including simultaneous data from Swift in the UV and soft X-ray band, from NuSTAR in the hard X-ray band and Fermi in the HE γ-ray band, together with contemporaneous data at very high energy (VHE) from H.E.S.S. and VERITAS, as well as archival data from WISE in the infrared. It should be noted that, although the VHE data are not strictly simultaneous, extreme-TeV blazars are known for their lack of variability in the VHE range. Some level of variability has so far only been observed for 1ES 0229+200 at a timescale of ∼10⁷ s and 1ES 1218+304 at a timescale of ∼10⁵ s.

The set of SEDs shows a clear continuity of the spectral shape between the HE and VHE range, except for the source 1ES 0414+009, which exhibits a highly unusual spectral upturn between the Fermi-LAT and H.E.S.S. spectra. Such a concave shape, if confirmed, may indicate the presence of an additional spectral component, beyond the one-zone description we are providing here. It is however plausible, as assumed by Costamante et al. (2018), that a temporal variation in the spectral shape in the HE band is responsible for this apparent upturn between the non-simultaneous datasets. Given these uncertainties, this source has thus been excluded from our study, reducing the sample to five sources.

As a first step, the values of the model parameters were determined using the analytical equations given in Sect. 2.5. A high Doppler factor of δ = 50 was chosen initially following the results of Costamante et al. (2018) and adjusted to lower or higher values if a good representation of the data could not be achieved. The estimated parameter values were used as a starting point for modelling with a numerical code that treats the radiative transfer and all relevant leptonic and hadronic emission processes (Cerruti et al. 2015; Zech et al. 2017). The contribution from the host galaxy was reproduced with a template from Silva et al. (1998), where the normalisation was adjusted to fit the infrared and optical data points. We have verified, analytically and with the time-dependent code by Dmytriiev et al. (2021), that in the resulting solutions, radiative cooling of the emitting particle population can be safely neglected, given the low magnetisation that is inherent to the shock acceleration scenario.

The co-acceleration model was declined into three scenarios to find the physically most meaningful representations of the selected SEDs. In a first scenario, we study a generic version with a large value of γ_e, min, mimicking electron preheating at a shock of large Lorentz factor, or in some unspecified primary dissipation region. In a second scenario, we study the case of a lower value of γ_e, min, as discussed earlier in the context of preheating at mildly relativistic internal or recollimation shocks. Finally, in a third scenario, we investigate the possibility of reacceleration at multiple shocks, focussing on the case of a single reacceleration at a secondary shock, sufficient to provide a good representation for the datasets under study, as will be seen.

In all cases, synchrotron radiation from the protons peaked at radio frequencies and was several orders of magnitude lower than the emission from electrons. It can thus be completely neglected when studying the SEDs. Only the kinetic energy of the protons needs to be considered when discussing the energetics of the source.

3.1. Scenario I: Generic shock model with large γ_e, min

We thus consider here a model with a large value of γ_e, min. In light of our discussion in Sect. 2, such a large value can either represent the consequence of electron preheating at a shock of large Lorentz factor or some pre-energisation in some dissipation process located well upstream of the shock where electron acceleration takes place. In the latter two-zone model, protons do not play any role at the shock; it would therefore apply equally well for a pure or high-multiplicity electron-positron composition.

In the former, a large value of the shock Lorentz factor (e.g., ∼δ), means that the jet encounters an obstacle moving at small velocity in the source rest frame (see Sect. 2.4.1). However, this requires the jet ram pressure to exceed that of the obstacle in the source frame, which is precisely opposite to the condition that allows the obstacle to penetrate the jet (e.g., Araudo et al. 2010; Barkov et al. 2010; Bosch-Ramon et al. 2012). Consequently, none of the above possibilities currently receives clear justification. For this reason, we regard this first scenario as a case study. To fix the parameters, we assume that γ_e, min is set by the shock Lorentz factor, in accordance with Sect. 2.3, and that this shock Lorentz factor is set equal to the Doppler factor δ, then slightly adjusted by constraints from the slope of the Fermi-LAT spectrum where necessary. This turned out to be necessary for one source only, 1ES 1218+304 (see below).

This choice leads to very high values of γ_e, min, of a few 10⁴, in agreement with previous modelling attempts that can be found in the literature, where γ_e, min is treated as a free parameter. We further assumed that the magnetisation in the emission region is the same as in the region upstream of the shock (i.e. σ_rad = σ). This choice of constraints removes the usual degeneracy of the SSC model. Ideally, following the procedure outlined in Sect. 2.5, only the Doppler factor δ is varied to adjust the model to the SED.

Without any further adjustments, this scenario leads to a satisfactory representation of the SED of 1ES 0229+200 for δ = 50 (cf. Fig. 1). Equally good solutions can be found for RGB J0710+591, 1ES 0347−121, and 1ES 1101−232 (cf. Figs. D.1–D.3) if, in addition to adjusting δ to values between 30 and 50, one also allows for an adjustment of the size of the emission region R, typically by a factor of ∼1.5. In the case of 1ES 1218+304, the Fermi-LAT spectrum is flatter than for the other sources, requiring a decrease in γ_e, min by an order of magnitude, in addition to adjustments of δ and R (cf. Fig. 2).

Fig. 1. SSC model for the SED of 1ES 0229+200 in scenario I. The data points at VHE are de-absorbed on the extragalactic background light (EBL) using the model by Franceschini et al. (2008) in this figure and in all following figures showing SEDs.

Fig. 2. SSC model for the SED of 1ES 1218+304 in scenario I.

Allowing for these adjustments, satisfactory solutions can be found for all five SEDs. The full application of the constraints of the shock co-acceleration model leads to unique solutions in this generic scenario. As can be seen from Table 1, the parameter values are very similar for the different sources: the Doppler factor of the emission region is large (up to δ = 60); the emission region has a size of a few 10¹⁶ cm, a standard value in SSC models for high-frequency-peaked BL Lac objects (HBLs); the magnetic field strength, in the milligauss range, and the magnetisation are very small, as has been found before for these extreme sources. In our framework, a low magnetisation σ ≪ 10⁻² is of course required to guarantee efficient acceleration.

In certain cases, for example the model of 1ES 1218+304 shown in Fig. 2, a mismatch between the model and the SED appears in the ultraviolet and infrared ranges. This may be explained with an additional synchrotron emission component from the extended jet, which is not represented in our models.

3.2. Scenario II: Co-acceleration on mildly relativistic shocks

In this second scenario, we consider the more physically grounded case of electron preheating and acceleration at internal (Sect. 2.4.1) or recollimation shocks (Sect. 2.4.2). In both cases, far smaller values of the shock Lorentz factor in the shock normal frame Γ_j < |n should be assumed than those applied in the generic scenario I.

For acceleration on an internal shock caused by a plasma blob with Lorentz factor Γ_b traversing the jet, the resulting shock Lorentz factor γ_sh ≡ Γ_j < |n is given by Eq. (4). Assuming a Doppler factor of δ for the observed emission from the blob, a mildly relativistic shock with Γ_j < |n ≃ 3 would result for a jet Lorentz factor Γ_j ≃ δ/12, when assuming a closely aligned jet, corresponding to Γ_b ≃ δ/2, such that Γ_j would be of the order of a few. Such a combination of a highly relativistic emission component responsible for HE emission and a mildly relativistic jet is frequently proposed to explain the discrepancy in observed bulk Lorentz factors in the radio and HE bands (Ghisellini et al. 2005; Henri & Pelletier 1991; Sol et al. 1989). In this case, shock Lorentz factors much higher than a few can only be achieved for unrealistically high blob Lorentz factors or for non-relativistic jets.

Considering the second case of acceleration on recollimation shocks, since the bulk Lorentz factor of the downstream, radiating plasma is Γ_j > |s ≥ δ/2, the shock Lorentz factor must verify Γ_j < |n ≃ Γ_j < |s/Γ_j > |s ≲ 2Γ_j < |s/δ (cf. Appendix B). If the bulk Lorentz factor of the plasma flow upstream of the recollimation shock Γ_j < |s is of a similar strength as the typical bulk Doppler factors derived from SED modelling, the shock Lorentz factor Γ_j < |n should be close to unity. The default value of Γ_j < |n = 3 that we impose in this second scenario for both cases to better accommodate the datasets is already relatively high, but it should still be acceptable in this case.

With this assumption, the minimum electron Lorentz factor is fixed to γ_e, min ≃ 600 Γ_j < |n ≃ 1800. The corresponding model fits imply different values of the physical parameters, following the procedure outlined in Sect. 2.5. Retaining bulk Doppler factors around δ ≈ 50, as for the solutions in scenario I, the decrease in γ_e, min for a fixed γ_e, max leads to a quadratic decrease in the inferred upstream magnetisation σ (cf. Eq. (3)). Since the particle density is inversely proportional to the magnetisation downstream of the shock (cf. Eq. (14)), the condition σ_rad = σ needs to be abandoned to avoid a strong increase in the particle density, which would lead to a Compton dominance in the SED that is not observed in extreme blazars.

An increase in the magnetisation in the region downstream of the shock by one or two orders of magnitude is not unexpected, due to the self-generation of magnetic turbulence at the shock. As discussed earlier, shock microphysics only imposes the condition σ ≤ σ_rad ≲ 0.01, with 0.01 representing the typical effective magnetisation at the shock, σ_rad the magnetisation in the radiation region and σ the pre-shock magnetisation. We also stress that the values of σ that are provided in Table 1 are only indicative, in the sense that they are derived from the respective values of γ_e, min and γ_e, max under the assumption that the maximal energy is set by the scattering constraint. As discussed earlier, this constraint applies formally to super-luminal shocks; it would be relaxed in sub-luminal shocks and, in such cases, the actual pre-shock magnetisation could be as large as σ_rad. Consequently, the true physical magnetisation to be regarded is σ_rad, of the order of 10⁻⁴–10⁻³ for the various sources in this model.

As can be seen from Figs. 3, D.4 and D.5, satisfactory representations of the SEDs can be found within this shock scenario for 1ES 1218+304, RGB J0710+591 and 1ES 0347−121. While for the first two, the representation of the SED is of comparable quality to the one in scenario I, for 1ES 0347−121 the Fermi-LAT flux is over-predicted for the lowest flux point, but given the uncertainties in the Fermi points, this solution can still be regarded as acceptable. It should be noted here that, as far as the synchrotron and SSC peak flux levels are concerned, the residual freedom in σ_rad introduces a degeneracy in the sense that similar solutions can be found by increasing δ and decreasing σ_rad or vice versa.

Fig. 3. SSC model for the SED of 1ES 1218+304 in scenario II.

In the case of 1ES 0229+200, the SSC component of the model is not sufficiently steep to pass through the γ-ray data points, as a result of the smaller value of γ_e, min. The estimated SSC peak frequency, which is one of the input parameters for the modelling, needs to be increased significantly to be able to approach the VHE data points in 1ES 0229+200. This leads to an increased radius of the emission region, of the order of a parsec, and a required jet power close to the Eddington limit. The magnetic field strength of this solution would be quite small, of the order of 10⁻⁵ G.

The final source in the sample, 1ES 1101−232, cannot be satisfactorily described either given the model constraints. In addition to a poor representation of the Fermi-LAT points, the infrared flux is overproduced by the model. If γ_e, min remains fixed at its default value, the only way around this problem is to increase δ to a value well above 100, in order to shift the low-energy break to sufficiently high frequencies.

Otherwise, satisfactory solutions for 1ES 0229+200 and 1ES 1101−232 can be obtained only if the shock Lorentz factor Γ_j < |n is increased to large values, thus recovering a scenario similar to I. These solutions are shown in Figs. D.6 and D.7 and the corresponding parameters are included in Table 1. To reproduce the measured VHE flux in 1ES 0229+200, a minimum value of Γ_j < |n ≈ 20 is necessary, while in the case of 1ES 1101−232, Γ_j < |n ≈ 10 leads to a sufficiently steep rise of the SED at low energies to respect the constraints from the optical/IR data.

Given these solutions, if acceleration were to take place at a mildly relativistic recollimation shock, this would lead to an initial bulk Lorentz factor of the jet of Γ_j < |s ≈ 10³ for 1ES 0229+200, while for 1ES 1101−232, it would result in Γ_j < |s ≈ 500. Such large jet Lorentz factors are problematic given that the usually assumed values for less extreme blazars, based on observations with very-long-baseline interferometry (VLBI) and on models of relativistic jets, are of the order of a few to a few tens (Pushkarev et al. 2009; Saikia et al. 2016).

For acceleration on internal shocks caused by blob-jet interaction, the required values of Γ_j < |n can be accommodated with jet Lorentz factors of a few, if the jet is not too closely aligned with the line of sight (i.e. for δ ≲ Γ_b).

3.3. Scenario III: Re-acceleration on a second shock

Despite the difficulties the recollimation shock scenario has in accounting for the SEDs with the hardest γ-ray spectra, it remains attractive as a natural explanation for continuous non-variable emission. For this reason, we explore here the possibility that particles are accelerated at multiple standing shocks. As will be seen, such a possibility may help understand the peculiar SEDs of the most emblematic extreme blazars 1ES 0229+200 and 1ES 1101−232.

More generally, the need for continuous acceleration or multiple acceleration episodes along the jet is strongly supported by the non-thermal emission from radio galaxies at large (kpc) distances from the core, habitually observed in the radio, optical and X-ray bands, and recently even in VHE γ-rays (H.E.S.S. Collaboration 2020). For simplicity, and to keep the model as economical as possible, we consider the simplest generalisation of scenario II, namely re-acceleration on a second shock of a particle distribution that had already been accelerated on an initial shock. As long as cooling remains negligible, re-acceleration leads to a hardening particle spectrum with a low-energy turnover at increasingly higher energies above the initial value of γ_e, min, as detailed in Sect. 2.4.3.

The effect on the SED can be roughly estimated by approximating the electron distribution after n re-accelerations with a hardening power law with increasing values of ${\gamma }_{\mathrm{e},\min }^{(n)}\simeq {\mathit{g}}^n{\gamma }_{\mathrm{e},\min }$. If the maximum energy is limited by the scattering constraint, it increases as well by gⁿ (i.e. ${\gamma }_{\mathrm{e},\max }^{(n)}\simeq {\mathit{g}}^n{\gamma }_{\mathrm{e},\max }$), provided the magnetisation remains comparable from shock to shock. It cannot, of course, exceed the synchrotron or age limit. Fits of a power law to the distribution given by Eq. (6), far from ${\gamma }_{\mathrm{e},\min }^{(n)}$ and assuming g = 2 and a power law with index s⁽⁰⁾ = 2.2 for the initial acceleration, yield indices s⁽¹⁾ ≃ 2.1 and s⁽²⁾ ≃ 2.0 for the first and second re-acceleration. In this approximation, the synchrotron and IC peak frequencies increase as ${\nu }_{\mathrm{s}}^{(n)}\propto {\mathit{g}}^{2n}{\nu }_{\mathrm{s}}$ and roughly ${\nu }_{\mathrm{IC}}^{(n)}\propto {\mathit{g}}^n{\nu }_{\mathrm{IC}}$, respectively.

The total emitted synchrotron power per unit volume and frequency at the frequency ν for the power-law electron distribution after n re-accelerations in a turbulent magnetic field is then given following Rybicki & Lightman (2004):

$$\begin{array}{cc}& P_{\mathrm{syn}}^{(n)}(\nu )\propto K_{\mathrm{e}}^{(n)}\frac{1}{s^{(n)}+1}\mathrm{\Gamma }(\frac{s^{(n)}}{4}+\frac{19}{12})\mathrm{\Gamma }(\frac{s^{(n)}}{4}-\frac{1}{12})\hfill \\ \hfill & \times {\left(\frac{2\pi \nu m_{\mathrm{e}}c}{3eB}\right)}^{-(s^{(n)}-1)/2},\hfill \end{array}$$(17)

with $K_{\mathrm{e}}^{(n)}$ the normalisation of the power-law distribution.

Keeping the electron number density and magnetic field strength constant, and estimating the peak frequency of the synchrotron component in the δ approximation as ${\nu }_s^{(n)}[\mathrm{Hz}]\simeq 3.7\times {10}^{6}B[\mathrm{G}]{{\gamma }_{\mathrm{e},\max }^{(n)}}^2$, the increase in the emitted peak power for the first re-acceleration is roughly

$$\begin{array}{c}\hfill \frac{{\nu }_s^{(1)}{P}_{\mathrm{syn}}^{(1)}({\nu }_s^{(1)})}{{\nu }_s^{(0)}{P}_{\mathrm{syn}}^{(0)}({\nu }_s^{(0)})}\simeq 7.1{\left(\frac{{\gamma }_{\mathrm{e},\max }}{{10}^6}\right)}^{0.1}.\end{array}$$(18)

The numerical factor in the above equation depends weakly on the value of γ_e, min, which we set again to 1.8 × 10³ here.

Under the same assumptions, the Compton dominance of the emission (i.e. the ratio of the luminosities of the IC and synchrotron components) increases with re-acceleration. Following Tavecchio et al. (1998), it can be written as

$$\begin{array}{cc}& {\kappa }^{(n)}\equiv {\left(\frac{L_{\mathrm{IC}}}{L_{\mathrm{syn}}}\right)}^{(n)}\propto f({\alpha }_1^{(n)},{\alpha }_2){\nu }_{\mathrm{s}}^{(n)}{P}_{\mathrm{syn}}^{(n)}({\nu }_{\mathrm{s}}^{(n)})\hfill \\ \hfill & \times {\left(\frac{3m_{\mathrm{e}}{c}^2}{4\times 3.7\times {10}^{6}h}\right)}^{(3-s^{(n)})/2}{\left(\frac{\delta }{B{{\gamma }_{\mathrm{e},\max }^{(n)}}^3}\right)}^{(3-s^{(n)})/2}.\hfill \end{array}$$(19)

The ratio of the Compton dominance between the first re-acceleration and the initial acceleration is then given by:

$$\begin{array}{c}\hfill \frac{{\kappa }^{(1)}}{{\kappa }^{(0)}}\simeq 2.6{\left(\frac{\delta }{50}\frac{1\mathrm{mG}}{B}\frac{{10}^6}{{\gamma }_{\mathrm{e},\max }}\right)}^{0.05}.\end{array}$$(20)

The actual spectral distribution of re-accelerated particles deviates from a pure power law, leading to an even more significant effect on the resulting SED. The evolution of a given electron population and its associated broadband emission from an initial shock acceleration to a first and second re-acceleration on similar shocks is shown in Figs. 4 and 5. For this simulation, it was assumed that the number of particles remains constant during the evolution and that radiative cooling can be neglected. The source parameters δ, R, and B were adapted from the solution for 1ES 0229+200 that will be presented below and were left constant for the consecutive accelerations. Under these admittedly simple assumptions, it can be seen that the resulting SED of the electron population is dominated by the last efficient shock acceleration.

Fig. 4. Particle spectra for the same electron and proton population after the initial shock acceleration and after a first and second re-acceleration on successive shocks.

Fig. 5. Observed SEDs due to SSC emission from the particle population shown in Fig. 4 after an initial shock acceleration and two re-accelerations. The SEDs correspond to a source at the redshift of 1ES 0229+200 with δ = 50.

Assuming that the emission from 1ES 0229+200 and 1ES 1101−232 is dominated by radiation of a particle population after a first re-acceleration, the SEDs of both sources can be well represented with parameters that are similar to scenario II; in particular a value of γ_e, min = 1.8 × 10³ is sufficient. To avoid an overestimate of the synchrotron peak frequency in the case of 1ES 1101−232, due to the harder particle spectrum, the input value ν_syn for the model was reduced by a factor of 2. The results of this scenario are shown in Figs. 6 and 7, with the parameter values given in Table 1. This third scenario also leads to a better representation of the SED of 1ES 0347−121 (cf. Fig. D.8).

Fig. 6. SSC model for the SED of 1ES 0229+200 in scenario III.

Fig. 7. SSC model for the SED of 1ES 1101−232 in scenario III.

4. Discussion

4.1. Implications for the nature and location of the emission region

The emission region described by the simple one-zone model may represent either a plasma blob moving with relativistic speed through the jet, into which accelerated particles are injected from a leading bow shock (similar to Kirk et al. 1998), or the region in the relativistic plasma flow downstream from a standing shock.

In both cases, one needs to explain why the emission is restricted mostly to a compact region, resembling the moving or standing VLBI radio knots observed in radio galaxies and blazars. Since in our scenarios synchrotron cooling is very slow, adiabatic cooling or particle escape into regions with even lower magnetic fields may be an explanation. As was discussed above, magnetisation is in fact expected to be amplified in the immediate neighbourhood of the shock. When particles escape from the emission region into the downstream jet, their emission may be shifted to lower fluxes and frequencies. This low-energy emission from the extended jet may be responsible for the radio flux, which we do not include in the SEDs shown here (cf. archival data in Costamante et al. 2018), and possibly contribute at infrared to ultraviolet frequencies, which may help to better reproduce the SED of 1ES 1218+304 in this range (cf. Fig. 3). Modelling of this extended emission component is beyond the scope of the present study.

The blob-in-jet scenario can accommodate reasonably well a baryon loaded medium with low magnetisation if the blob interacts with the jet at large distances from the jet base. However, the blob emission timescale is constrained both by the expansion of the blob and its slowdown, due to its interaction with the jet. We assume that the expansion of the blob takes place at a lateral velocity v_⊥ in the source frame, with v_⊥/c ∼ 𝒪(Θ_j) if the expansion follows that of the jet (e.g., Sbarrato et al. 2011). Then the available dynamical time in the blob frame is $t_{\text{dyn}}^{\prime }~R/(v_{\perp }{\mathrm{\Gamma }}_{\text{b}})$, corresponding to an observer timescale of $\mathrm{\Delta }{t}_{\text{obs}}=t_{\text{dyn}}^{\prime }/\delta ~R/\left(v_{\perp }{\mathrm{\Gamma }}_{\text{b}}\delta \right)$. If v_⊥/c ∼ Θ_j, the rather short dynamical timescale would imply substantial variations in the light curve on timescales shorter than a year. If the jet opening angle is smaller than the usual rough assumption (i.e. in the case of Θ_j ≪ 1/Γ as discussed below), one could get a result closer to the observed slow or seemingly absent variability. In a transverse structured jet, the observed opening angle of the outer jet may also not be directly constraining the size of a blob that is propagating in an inner ‘spine’.

The blob slows down in the jet frame once it has swept up a mass $M_{\text{j}}~E_{\text{b}}/({\mathrm{\Gamma }}_{\text{b}|\text{j}}^2{c}^2)$, where E_b represents the blob energy (in the jet frame), and Γ_b|j the relative blob-jet Lorentz factor. On a timescale $t_{\text{dec}}^{\prime }$ (blob frame), the swept-up mass is $M_{\text{j}}~\pi R^2{n}_{\text{j}}{m}_{\text{p}}{\mathrm{\Gamma }}_{\text{b}|\text{j}}c{t}_{\text{dec}}^{\prime }$, hence the deceleration time can be written $t_{\text{dec}}^{\prime }~\left[u_{\text{b}}/\left({\mathrm{\Gamma }}_{\text{b}|\text{j}}{n}_{\text{j}}{m}_{\text{p}}{c}^2\right)\right]R/c$, with u_b the blob proper energy density. Consequently, $t_{\text{dec}}^{\prime }$ is generally larger than $t_{\text{dyn}}^{\prime }$, because the blob ram pressure ${\mathrm{\Gamma }}_{\text{b}|\text{j}}^2{u}_{\text{b}}$ is assumed larger than n_jm_pc², and Γ_b|j is of the order of a few at most.

The alternative scenario of acceleration on standing shocks, arising for example from recollimation, may account more naturally for a continuous emission. In the standard picture of an initially Poynting-flux dominated jet, a first recollimation shock is assumed to occur at the transition between the magnetically dominated parabolic base of the jet and a kinetically dominated conical jet at the pc-scale (e.g., Potter 2018). Oblique shocks due to collimation of the relativistic jet by the pressure of the surrounding medium have been considered by, for example, Bromberg & Levinson (2009) to account for the difference between large Doppler factors needed to explain TeV emission and the much smaller values inferred from unification models and super-luminal motion of radio knots. In this framework, acceleration and emission up to VHEs can occur up to several 10 pc from the base of the jet, from very compact emission regions with extensions down to ∼10⁻³ times the distance scale, in agreement thus with the size scale of the emission regions in our models.

Direct observational evidence for the existence of such recollimation shocks has so far only been found from VLBI data in a few radio-loud AGN (e.g., Hada et al. 2018, and references therein), and most prominently in the most nearby radio galaxy M 87 at the ‘HST-1’ radio knot (Asada & Nakamura 2012; Nakamura & Asada 2013), which appears at a de-projected distance of approximately 300 pc from the core, with a radius of the order of 1 pc (Asada & Nakamura 2012).

The typical extensions of the emission regions in our models are of the order of 0.01 pc. To explain the HE emission from extreme-TeV blazars, much larger regions can be ruled out due to the energy requirements, as discussed in the case of 1ES 0229+200 for our scenario II. Besides the above-mentioned predictions from Bromberg & Levinson (2009), an order-of-magnitude estimate of the location of the emission regions in our models can be made by comparing their radii with the assumed jet opening angle Θ_j, since Θ_j ∼ R/d, with d the distance of the emission region from the base of the jet. Typical values for Θ_j, based on VLBI data from the MOJAVE programme for γ-loud blazars are found to follow roughly Θ_j ≃ 0.3/Γ_{j|< s} (Pushkarev et al. 2009), while Zdziarski et al. (2015) find Θ_j ≃ 0.1/Γ_j in a different sample. The latter note however that the opening angle is related to the magnetisation as Θ_j ≃ sσ^1/2/Γ_{j|< s} with s ≲ 1, which would suggest a relation closer to Θ_j ≃ 0.01/Γ_{j|< s} for the typical values of σ found in our scenarios². These estimates lead to a range of distances from the core of the order of one to several tens of parsecs.

Re-acceleration on successive standing shocks requires the distance between two shocks to be small compared to the cooling time of particles inside the plasma flow. Following Eq. (2), the most energetic electrons in our models will not significantly cool at distance scales up to about 0.1 Γ_{j|< s} pc (i.e. a few parsecs). This estimate is based on the magnetic field strength in the emission region, which is expected to be amplified compared to regions far from the recollimation shocks; the actual cooling distance may thus be far greater. Formation of internal shocks due to reflections of the converging recollimation shock, proposed by Bromberg & Levinson (2009), may constitute alternative sites of re-acceleration.

4.2. Implications for jet power, jet content, and magnetisation

Contrary to other lepto-hadronic models in the literature, which try to explain the observed SEDs as a combination of radiative emission from electrons and protons, emission from protons is negligible in the scenarios proposed here, due to the low magnetic field strengths and small maximum proton Lorentz factors γ_p, max. Typical lepto-hadronic models achieve a significant radiative contribution from protons by assuming an acceleration process that is much more efficient for protons than for electrons, leading to u_p ≫ u_e and γ_p, max ≫ γ_e, max, combined with high magnetic field strengths, of the order of 1 G to a few 100 G in the case of proton-synchrotron dominated emission. In such models, the jet power is strongly dominated by the contribution from either the magnetic field or the relativistic protons, depending on the dominant hadronic emission mechanism, and is generally close to or even above the Eddington limit (e.g., Mücke & Protheroe 2001; Petropoulou & Mastichiadis 2012; Cerruti et al. 2015, ...).

In our solutions, the magnetic field strength is several orders of magnitude smaller, the number densities of accelerated electrons and protons are equal (n_e = n_p), and the jet power is typically a factor 10² to 10³ below the Eddington limit, closer to the results for BL Lac objects seen in the leptonic framework (e.g., Ghisellini et al. 2011). By design, emission emerges from a region of the jet with a low magnetisation σ ≪ 10⁻².

Studies of luminous, FSRQ-type blazars show that, at least for those objects, the assumption of a pure electron-proton jet is problematic, since required jet powers are very large compared to the available accretion power and compared to the power that is injected into radio lobes. The addition of an electron-positron pair plasma can lower the required jet power (e.g., Sikora 2016). If this is the case, n_e > n_p, which would change our constraint on the normalisation of the electron spectrum K_e (Eq. (15)). The value of γ_e, min is expected to decrease by the factor of electron multiplicity. A large multiplicity factor would thus exclude a good representation for the scenarios without substantial re-acceleration.

The magnitude of the magnetisation σ of blazar jets is still an open question. When assuming the Blandford-Znajek process for jet launching, this implies an initially strongly Poynting-flux dominated jet, in which magnetic energy is converted to kinetic energy over subparsec to parsec scales. During this conversion, the bulk Lorentz factor increases and the magnetisation drops to unity and below. In the standard picture, this conversion may proceed through differential collimation down to equipartition (i.e. σ ≃ 1), where it becomes inefficient.

Further conversion may be driven by non-standard mechanisms, such as reconnection on magnetohydrodynamic (MHD) instabilities (Komissarov 2011). In the latter case, a magnetisation σ ≪ 1 seems reachable even at small distances from the base of the jet, where most of the radiation is expected to be produced.

In this case, Zdziarski et al. (2015) show that blazar jets can be powered by the Blandford-Znajek mechanism, originate from magnetically arrested disks and can still have magnetisations σ ≪ 1 at the radio core. They also find that in this framework jets are relatively heavy, with their mass dominated by baryons.

Such reconnection processes could lead to particle acceleration. However, if the magnetic field strength in the transition region around σ ≲ 1 does not exceed values of the order of 1 G, the emission from those particles would contribute mainly to the extended radio and optical emission of the jet, not to the high-frequency part, because the average electron Lorentz factor decreases with decreasing σ and the produced spectra become quite soft for σ ≲ 1 (see Appendix A). Such an additional emission component might actually help explain the discrepancy between our model and certain SEDs at low energies.

A very rough comparison between the synchrotron peak emission close to such a transition region and the one from the HE emission zone can be made by considering two spherical regions with, in one case, B_trans ≲ 1 G, σ_trans ≲ 1 and a steep particle spectrum with γ_e, min ∼ 10² and, in the other case, average conditions from our results from scenario II. Estimating the peak luminosities using Eq. (17) leads to a ratio of roughly L_syn, trans/L_syn, he ∼ 50 (R_trans/10¹⁶ cm)³(B_trans/1 G)⁴(1/σ_trans), with R_trans the radius of the transition region. For sufficiently small values of B_trans and R_trans, this emission may be hidden below that of the host galaxy and HE region.

The luminosity of the IC peak from such a transition region would be much weaker than the synchrotron luminosity due to the steep particle spectrum. A more thorough treatment of this question would require simulations with a multi-zone jet emission model including a realistic description of the evolving reconnection regions, which is well beyond the scope of the current work.

In alternative scenarios where jet launching is described through the Blandford-Payne mechanism, a lower initial magnetisation of the jet may arise more naturally (e.g., Tzeferacos et al. 2009, for simulations in the non-relativistic limit).

Direct measurements of the magnetic field strength B in relativistic jets of blazars and radio galaxies are difficult. They usually rely on measures of the core shift, that is to say, the frequency-dependent position along the jet where synchrotron self-absorption becomes negligible and emission from the VLBI radio core can first be observed. To determine the value of B, usually equipartition between magnetic and particle energy densities is assumed and values of the order of 100 mG are found at a distance of 1 pc from the base of the jet (O’Sullivan & Gabuzda 2010; Kutkin et al. 2014), while Agarwal et al. (2017) find values of the order of 1 mG at the more distant radio cores. Uncertainties and fluctuations in these measurements are however large and deviation from the assumption of equipartition will lead to different values.

4.3. Applicability to other blazar types

The SEDs of less extreme HBLs are well represented with smaller values of γ_e, min, typically of the order of 1–100. These more common objects present at least an order of magnitude lower synchrotron and SSC peak frequencies than the extreme-TeV blazars, which translates into lower values of γ_e, max and δ, while the magnetic field strength is usually found to be at least an order of magnitude higher.

In the present framework, the continuous emission from such objects could in principle be accounted for with lower Lorentz factors in the shock normal frame. A higher electron multiplicity would also lead to lower values of γ_e, min. Lower values of δ may be due to smaller Lorentz factors in the post-shock region or larger viewing angles of these objects. To reproduce the wider synchrotron and SSC bumps seen in many HBLs, particle spectra require in general a slightly more complex description, usually including a spectral break that cannot be ascribed to radiative cooling of a single electron population alone. These features could possibly be reproduced with a combination of emission from several standing shocks (e.g., a strong first shock) followed by a weaker second shock. In the blob-in-jet scenario, the assumption of a uniform particle density and magnetic field strength in a spherical blob may be too approximate for a realistic description of the particle spectra, hence requiring additional degrees of freedom in the form of an ad hoc spectral break. The acceleration process may also be more complex than what is assumed in our current model.

Contrary to extreme TeV blazars, flux and spectral variability are ubiquitous in HBLs, at timescales that can be as short as a few minutes in the observer frame at TeV energies. Additional mechanisms, possibly the interaction of moving shocks with successive standing shocks (e.g., Hervet et al. 2019; Fichet de Clairfontaine et al. 2021, and references therein) need to be considered to account for the appearance of high flux states and flaring events in blazars and radio galaxies. Alternatively, depending on the magnetisation in the emission region, scenarios based on particle acceleration through reconnection or turbulence may provide a suitable description for the variable emission from such objects.

5. Conclusions

Ascribing the non-thermal multi-wavelength emission from extreme-TeV blazars to the acceleration of a proton-electron plasma on mildly relativistic shocks introduces additional constraints, arising from the microphysics in the shock region, on the standard one-zone SSC model that is generally applied to describe SEDs of BL Lac objects. A more physically motivated description of the population of radiating electrons can thus be achieved. The co-acceleration of electrons and protons results in a transfer of energy to the electrons, thus shifting their spectral distribution to higher energies. This effect provides a natural explanation for the high values of the minimum electron Lorentz factors required to model the SEDs of extreme-TeV blazars.

A first generic model, allowing for a very high bulk Lorentz factor of the jet plasma upstream of the shock, leads to unique solutions for the set of five extreme-TeV blazar SEDs under study. In a second, more realistic scenario, exploring the two cases of acceleration on recollimation shocks or on the shock wave caused by a plasma blob propagating inside the jet, a good representation is found for most of the SEDs. High minimum Lorentz factors of the electron distribution of γ_e, min ≳ 10³ arise as a natural consequence of the assumed microphysics.

For the SEDs of the two objects with the hardest γ-ray spectra, 1ES 0229+200 and 1ES 1101−232, good representations are found for acceleration on moving shocks, but also for the acceleration on recollimation shocks when assuming that the particle distribution that dominates the emission results from re-acceleration on a second shock. In the latter case, the spectral hardening and the increase in the minimum and maximum Lorentz factors arising from the re-acceleration of an initial power-law electron distribution on a second shock can satisfactorily explain the extreme values of the synchrotron and SSC peak frequencies and the hard γ-ray spectrum.

Acknowledgments

The authors wish to thank Luigi Costamante for providing the data points of the SEDs and Zakaria Meliani for useful exchanges.

References

Appendix A: Electron heating in reconnection layers or magnetised turbulence

To produce electron spectra with large ${\overline{\gamma }}_e$, reconnection or turbulent acceleration must proceed in the relativistic regime, because they need to dissipate an energy reservoir $\sim {\overline{\gamma }}_e/{\gamma }_0$ times larger than the initial electron energy content, as discussed earlier. In a pure pair plasma, this would mean an initial magnetisation $\sigma \sim {\overline{\gamma }}_e/{\gamma }_0\gg 1$, clearly at odds with the apparent requirement of a sub-equipartition magnetic field in the radiation region. As mentioned earlier, one may contemplate the possibility that such values are reached at the core of the jet in a region of high magnetisation, provided the radiation produced in this region does not dominate. Such dissipation generally leads, however, to an approximate equipartition between particles and magnetic fields, leading to a magnetisation σ ∼ 𝒪(1), which remains too large to account for the spectral properties of extreme blazars, whence the need for some additional mechanism to bring down this magnetisation to the desired sub-equipartition level by the time the electron power law is shaped.

In an electron-proton plasma, the initial energy density is carried by the protons, hence reconnection can proceed in the relativistic regime from the point of view of electrons (σ_e > 1, but in the sub-relativistic regime from the point of view of ions (σ_p ≲ 1). Here, σ_e (respectively, σ_p) denotes the electron (respectively, ion) magnetisation, that is, the ratio of the energy density in the magnetic field to the electron (ion) initial energy density. The total magnetisation σ is related to these two quantities via ${\sigma }^{-1}={\sigma }_e^{-1}+{\sigma }_p^{-1}$. The simulations of Petropoulou et al. (2019), although conducted in the regime σ ≳ 1, suggest ${\overline{\gamma }}_e\sim 1+0.03{\sigma }^{3/2}{m}_p/m_e$ for an initially cold electron population, hence values σ ≳ 10 appear needed to reach ${\overline{\gamma }}_e\sim {10}^3$. This result remains essentially unchanged if the electrons are initially hot (T_e ≫ m_e). At lower magnetisations, the electron spectra take a more complex form, composed by part of their initial Maxwellian distribution at low Lorentz factors and a non-thermal population emerging out of a secondary Maxwellian at large Lorentz factor ${\gamma }_e^{\text{hi}}$. From the results of Ball et al. (2018) and Werner et al. (2018), we infer ${\gamma }_e^{\text{hi}}~\left(0.1-0.3\right)\sigma m_p/m_e$. Since ${\overline{\gamma }}_e\lesssim {\gamma }_e^{\mathrm{hi}}$ for such electron spectra, we recover the above result: in other words, σ ≫ 1 is required to obtain ${\overline{\gamma }}_e\sim {10}^3$.

Werner et al. (2018) measure the fraction of dissipated energy q_e = e_e/(e_e + e_p), where e_e and e_p are understood as kinetic energies. They obtain q_e ∼ 0.25 + 0.25σ^1/2/(10+σ)^1/2 in the regime σ > 0.01, which, when combined with ${\overline{\gamma }}_p\simeq 1+0.5\sigma$ at σ ≲ 1 (Ball et al. 2018), gives ${\overline{\gamma }}_e\simeq 0.1\sigma m_p/m_e$ and thereby confirms the previous estimate.

Finally, we note that the slope of the accelerated spectrum depends rather sensitively on the magnetisation, for example s ≃ 4 for σ ≲ 1 and becoming as hard as s ∼ 2 for σ ≫ 1. All this clearly indicates that reconnection does not fulfil the requirements for a successful model for extreme blazars under our present assumptions.

The partition of energy between electrons and ions in MHD turbulence is a rather convoluted problem, which depends on a number of parameters, not only the initial magnetisation and species temperatures, but also the nature (compressive or Alfvénic) of the turbulence itself (e.g. Kawazura et al. 2020). As of today, there exist only a few kinetic simulations of ion-electron turbulence in the regime of interest and we rely on the results of Zhdankin et al. (2019). These authors have performed a scan of simulations which span a broad range of magnetisations, from σ_p ∼ 10⁻³ (hence, σ_e ∼ 1) up to the fully relativistic regime σ_p ≳ 1 (σ_e ∼ 10³). They observe that ions tend to be preferentially heated over electrons and that the amount of electron heating is found to scale with the initial temperature. For an initial electron temperature T_e ≃ m_e, the electrons acquire about 10% of the ion energy gain. If γ_p − 1 ∼ σ, as expected for near-complete transfer of the magnetic energy into particles, this means ${\overline{\gamma }}_e\simeq 0.1\sigma m_p/m_e$, a result comparable to what is observed in reconnection layers. Hence, ${\overline{\gamma }}_e\gtrsim {10}^3$ can be achieved for large magnetisations only.

Appendix B: Jump conditions at a recollimation shock

Here we concentrate on the shock crossing conditions and rewrite them in a simple way, in order to derive the post-shock flow velocity and to keep the discussion self-contained. We do so, in particular, by de-boosting to the so-called ‘shock normal frame’ ℛ_n (Begelman & Kirk 1990) in which the flow moves along the shock normal. This means de-boosting by the component of the velocity which is transverse to this shock normal (equivalently, parallel to the shock surface).

We assume that the surface of the recollimation shock makes an angle α with respect to the jet axis. The flow direction, at angle θ_j< with respect to this jet axis, thus forms an angle π/2 − (θ_j<−α) with respect to the shock normal. The symbol _< is used to distinguish pre-shock from post-shock quantities (indexed with _>). If the shock front is stationary in the source rest frame ℛ_s, then the shock normal frame ℛ_n moves at velocity β_n|s with respect to ℛ_s, where

$$\begin{array}{c}\hfill {\mathit{\beta }}_{\mathbf{n}\mathbf{|}\mathbf{s}}={\mathit{\beta }}_{\mathbf{j}\mathbf{<}}-({\mathit{\beta }}_{\mathbf{j}\mathbf{<}}·{\mathit{n}}_{\mathbf{|}\mathbf{s}}){\mathit{n}}_{\mathbf{|}\mathbf{s}}.\end{array}$$(B.1)

All throughout n_|s represents the (unit) direction of the normal to the shock front in the source rest frame. This gives β_n|s = β_j<cos(θ_j<−α), hence ${\mathrm{\Gamma }}_{\text{n}|\text{s}}={\left[1-{\beta }_{\text{j}<}^2{\cos }^2\left({\theta }_{\text{j}<}-\alpha \right)\right]}^{-1/2}$.

To derive the effective velocity of the jet in this shock normal frame, we note that the component of the jet four-velocity along the shock normal is preserved by the de-boost, hence u_j < |n = Γ_j<β_j<sin(θ_j<−α). The relative Lorentz factor between the recollimation shock and the incoming plasma in the shock normal frame, which matters for acceleration purposes, can thus be written

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{j}<|\mathrm{n}}=\sqrt{1+{\mathrm{\Gamma }}_{\mathrm{j}<}^2{\beta }_{\mathrm{j}<}^2{sin}^2({\theta }_{\mathrm{j}<}-\alpha )}=\frac{{\mathrm{\Gamma }}_{\mathrm{j}<}}{{\mathrm{\Gamma }}_{\mathrm{n}|\mathrm{s}}}.\end{array}$$(B.2)

If Γ_j<|θ_j<−α| ≪ 1, the recollimation shock is sub-relativistic in the ℛ_n frame, otherwise it is relativistic. The post-shock downstream velocity can be approximated (in ℛ_n) as β_j > |n ≃ β_j < n/r in the weakly magnetised limit, with r the compression ratio. Hereafter we assume a mildly relativistic shock with r ≃ 3, implying β_j > |n ≃ 0.3, hence Γ_j > |n ∼ 1. For reference, the exact shock jump conditions for a strong hydrodynamical shock give β_j > |n = 0.29 if Γ_j < |n = 3, β_j > |n = 0.26 if Γ_j < |n = 2.

Consequently, boosting back to the source frame gives a post-shock plasma four-velocity

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathbf{j}\mathbf{>}\mathbf{|}\mathbf{s}}={\mathrm{\Gamma }}_{\mathrm{j}>|\mathrm{n}}{\beta }_{\mathrm{j}>|\mathrm{n}}{\mathit{n}}_{\mathbf{|}\mathbf{s}}+{\mathrm{\Gamma }}_{\mathrm{j}>|\mathrm{n}}{\mathrm{\Gamma }}_{\mathrm{n}|\mathrm{s}}{\mathit{\beta }}_{\mathbf{n}\mathbf{|}\mathbf{s}},\end{array}$$(B.3)

with corresponding Lorentz factor

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{j}>|\mathrm{s}}={\mathrm{\Gamma }}_{\mathrm{j}>|\mathrm{n}}{\mathrm{\Gamma }}_{\mathrm{n}|\mathrm{s}}\simeq \frac{{\mathrm{\Gamma }}_{\mathrm{j}<}}{\sqrt{1+{\beta }_{\mathrm{j}<}^2{\mathrm{\Gamma }}_{\mathrm{j}<}{sin}^2({\theta }_{\mathrm{j}<}-\alpha )}}.\end{array}$$(B.4)

If Γ_j<|θ_j<−α| ≪ 1, the shock is weak, as noted above, and the post-shock plasma flows nearly along the shock surface. In the opposite limit, the recollimation shock is ultra-relativistic in ℛ_n (see Eq. (B.2)), the shock is strong and the post-shock flows at sub- or mildly relativistic speeds in ℛ_s (see Eq. (B.4)). To describe the intermediate limit, write |θ_j<−α| = κ/Γ_j< with κ ∼ 𝒪(1). Then ${\mathrm{\Gamma }}_{\mathrm{j}<|\mathrm{n}}\simeq \sqrt{1+{\kappa }^2}$ and ${\mathrm{\Gamma }}_{\mathrm{j}>|\mathrm{s}}\simeq {\mathrm{\Gamma }}_{\mathrm{j}<}/\sqrt{1+{\kappa }^2}$.

We also note that Eqs. (B.2) and (B.4) give Γ_j > |sΓ_j < |n  =  Γ_j < |sΓ_j > |n  ≃  Γ_j < |s. From the point of view of modelling, we recall that the post-shock Lorentz factor in the source rest frame (Γ_j > |s) controls the boost factor from the emission region, while the shock Lorentz factor in the shock normal frame (Γ_j < |n) controls the amount of electron heating (i.e. γ_{e, min} ≃ 600Γ_j<|n).

As the flow passes through the recollimation shock, it is refracted by an angle θ_j> − θ_j< (in the source rest frame), which can be obtained from the post-shock and pre-shock velocities. In the small-angle approximation, meaning |θ_j<|, |θ_j>|, |α|≪1, we obtain

$$\begin{array}{c}\hfill |{\theta }_{\mathrm{j}>}-{\theta }_{\mathrm{j}<}|\simeq \frac{|{\theta }_{\mathrm{j}<}-\alpha |}{\sqrt{3}}.\end{array}$$(B.5)

Therefore, if |θ_j< − α| ∼ 𝒪(1/Γ_j<), the deflection through the recollimation can lead to a substantial change in the Doppler beaming from pre-shock to post-shock values. We ignore this effect in our models, however.

Appendix C: Steady emission at a recollimation shock

As explained in Sikora et al. (1997), the Doppler amplification of the co-moving luminosity takes different forms depending on the geometric situation at hand: a moving blob of matter, or a steady jet. Both are reconciled once the steady jet is interpreted as a train of blobs moving one next to the other. The difference in amplification formula is then related to the actual number of blobs that is seen by an observer at a given time, as a result of relativistic motion.

For one moving blob, the general formula is:

$$\begin{array}{c}\hfill F_{\nu }=\frac{1+z}{d_{\mathrm{L}}^2}{\frac{\mathrm{d}E}{\mathrm{d}{t}_{\mathrm{obs}}\mathrm{d}\nu \mathrm{d}\mathrm{\Omega }}|}_{|\mathrm{s}}=\frac{1+z}{d_{\mathrm{L}}^2}{\delta }^3{\displaystyle \int }\mathrm{d}{V}^{\prime }{j}_{{\nu }^{\prime }}^{\prime }.\end{array}$$(C.1)

The subscript _|s indicates that the quantities are evaluated in the source rest frame. Henceforth, this subscript is dropped, and quantities written in the co-moving blob or jet frame are primed.

In the case of a steady jet, the standard formula reads

$$\begin{array}{c}\hfill F_{\nu }=\frac{1+z}{d_{\mathrm{L}}^2}\frac{{\delta }^2}{{\mathrm{\Gamma }}_{\mathrm{b}}}{\displaystyle \int }\mathrm{d}{V}^{\prime }{j}_{{\nu }^{\prime }}^{\prime }\end{array}$$(C.2)

because dt_obs = dt, versus dt_obs = dt(1−β_bμ) for the moving blob. However, describing the steady jet as a train of blobs leads to Sikora et al. (1997):

$$\begin{array}{c}\hfill F_{\nu }=\frac{1+z}{d_{\mathrm{L}}^2}{\delta }^3{N}_{\mathrm{b},\mathrm{eff}}{\displaystyle {\int }_{V_1^{\prime }}}\mathrm{d}{V}^{\prime }{j}_{{\nu }^{\prime }}^{\prime },\end{array}$$(C.3)

where $V_1^{\prime }$ denotes the co-moving volume of one blob, of co-moving size r′. The total volume of the emission region (in the source frame) is assumed to be composed of N_b blobs at a given time. Because the emission time of a given blob is compressed by (1 − β_bμ) = 1/(Γ_bδ) in the observer direction, the effective number of blobs that the observer sees at any given time is not N_b, but only N_b,  eff = N_b/(Γ_bδ). Then both descriptions indeed match one another, given that V = V′/Γ_b.

The total number N_b is fixed by the length scale over which emission lasts (in the source rest frame) for a given blob, that is, d_∥ (the length of the train of blobs in the source rest frame). The co-moving duration over which a blob is active is δt′=d_∥/(Γ_bβ_b) because δt′=δt/Γ_b and d_∥ = β_bδt (δt denoting the emission time in the source rest frame). The apparent width of a blob is δd_∥ = r′/Γ_b, hence

$$\begin{array}{c}\hfill N_{\mathrm{b}}=\frac{d_{\Vert }}{\delta d_{\Vert }}\simeq {\mathrm{\Gamma }}_{\mathrm{b}}^2{\beta }_{\mathrm{b}}\frac{\delta t^{\prime }}{r^{\prime }},\end{array}$$(C.4)

and consequently N_b,  eff ∼ δt′/r′. If δt′∼r′, as in the usual blob model, then N_b,  eff ∼ 1, hence the usual blob amplification formula can be used to describe the emission of the steady jet. The light curve will be roughly constant, because once a given blob has gone out of the (time) field of view of the observer, another has set in.

In principle, d_∥ can be fixed by other considerations than the assumption δt′∼r′. Consider for instance the following situation, which may apply to the case of a standing shock in the flow. One may consider that the magnetic field maintains a relatively constant amplitude in a region of extent d_∥, but then decays further on due to expansion, say as a some power law of d. If the magnetic energy density decays faster than d⁻¹, then it can be shown that particles do not actually cool in the decaying part, hence the emission is truly concentrated in the region in which the magnetic field strength is roughly constant. However, if r′ is set by the transverse size of the structure of the recollimation shock, d_∥ ∼ r′/Θ_j and Θ_j ∼ 1/Γ_b, one recovers N_b,  eff ∼ 1.

Finally, one cannot exclude altogether that N_b,  eff  ∼  a few, meaning that the observed flux would be the conjugation of several blobs at the same time. Nevertheless, adding the emission of several blobs would be similar to assuming an elongated emission region with a smaller radius, which would not fundamentally change the resulting emission.

Appendix D: Supplementary figures and parameters of the shock emission model

Solutions for the SEDs of the sources RGB J0710+591, 1ES 0347-121 and 1ES 1101-232 for the generic model (scenario I) are shown in Figs. D.1 to D.3. For the same sources, solutions of the recollimation-shock model (scenario II) are shown in Figs. D.4 to D.6. Figures D.8 and 7 show the solutions with a particle population that has been re-accelerated on a second shock (scenario III) for the sources 1ES 0347-121 and 1ES 1101-232, which feature very hard spectra in the Fermi-LAT range.

Fig. D.1. SSC model for the SED of RGB J0710+591 in scenario I.

Observables used as an input in the search for model solutions in all three scenarios, as well as parameters of the proton populations, are shown in Table D.1.

Fig. D.2. SSC model for the SED of 1ES 0347-121 in scenario I.

Fig. D.3. SSC model for the SED of 1ES 1101-232 in scenario I.

Fig. D.4. SSC model for the SED of RGB J0710+591 in scenario II.

Fig. D.5. SSC model for the SED of 1ES 0347-121 in scenario II.

Fig. D.6. SSC model for the SED of 1ES 1101-232 in scenario II.

Fig. D.7. SSC model for the SED of 1ES 0229+200 in scenario II.

Fig. D.8. SSC model for the SED of 1ES 0347-121 in scenario III.

All Tables

All Figures

	Fig. 1. SSC model for the SED of 1ES 0229+200 in scenario I. The data points at VHE are de-absorbed on the extragalactic background light (EBL) using the model by Franceschini et al. (2008) in this figure and in all following figures showing SEDs.
In the text

	Fig. 2. SSC model for the SED of 1ES 1218+304 in scenario I.
In the text

	Fig. 3. SSC model for the SED of 1ES 1218+304 in scenario II.
In the text

	Fig. 4. Particle spectra for the same electron and proton population after the initial shock acceleration and after a first and second re-acceleration on successive shocks.
In the text

	Fig. 5. Observed SEDs due to SSC emission from the particle population shown in Fig. 4 after an initial shock acceleration and two re-accelerations. The SEDs correspond to a source at the redshift of 1ES 0229+200 with δ = 50.
In the text

	Fig. 6. SSC model for the SED of 1ES 0229+200 in scenario III.
In the text

	Fig. 7. SSC model for the SED of 1ES 1101−232 in scenario III.
In the text

	Fig. D.1. SSC model for the SED of RGB J0710+591 in scenario I.
In the text

	Fig. D.2. SSC model for the SED of 1ES 0347-121 in scenario I.
In the text

	Fig. D.3. SSC model for the SED of 1ES 1101-232 in scenario I.
In the text

	Fig. D.4. SSC model for the SED of RGB J0710+591 in scenario II.
In the text

	Fig. D.5. SSC model for the SED of 1ES 0347-121 in scenario II.
In the text

	Fig. D.6. SSC model for the SED of 1ES 1101-232 in scenario II.
In the text

	Fig. D.7. SSC model for the SED of 1ES 0229+200 in scenario II.
In the text

	Fig. D.8. SSC model for the SED of 1ES 0347-121 in scenario III.
In the text