Issue 
A&A
Volume 654, October 2021



Article Number  A96  
Number of page(s)  18  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202141062  
Published online  19 October 2021 
Electronproton coacceleration on relativistic shocks in extremeTeV blazars
^{1}
Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France
email: andreas.zech@obspm.fr
^{2}
Institut d’Astrophysique de Paris, CNRS – Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France
email: lemoine@iap.fr
Received:
12
April
2021
Accepted:
31
July
2021
Aims. The multiwavelength emission from a newly identified population of ‘extremeTeV’ blazars, with Compton peak frequencies around 1 TeV, is difficult to interpret with standard onezone emission models. Large values of the minimum electron Lorentz factor and quite low magnetisation values seem to be required.
Methods. We propose a scenario where protons and electrons are coaccelerated on internal or recollimation shocks inside the relativistic jet. In this situation, energy is transferred from the protons to the electrons in the shock transition layer, leading naturally to a high minimum Lorentz factor for the latter. A low magnetisation favours the acceleration of particles in relativistic shocks.
Results. The shock coacceleration scenario provides additional constraints on the set of parameters of a standard onezone leptohadronic emission model, reducing its degeneracy. Values of the magnetic field strength of a few mG and minimum electron Lorentz factors of 10^{3} to 10^{4}, required to provide a satisfactory description of the observed spectral energy distributions of extreme blazars, result here from first principles. While acceleration on a single standing shock is sufficient to reproduce the emission of most of the extremeTeV sources we have examined, reacceleration on a second shock appears needed for those objects with the hardest γray spectra. Emission from the accelerated proton population, with the same number density as the electrons but in a lower range of Lorentz factors, is strongly suppressed. Satisfactory selfconsistent representations were found for the most prominent representatives of this new blazar class.
Key words: acceleration of particles / radiation mechanisms: nonthermal / BL Lacertae objects: individual: 1ES 0229+200 / BL Lacertae objects: individual: 1ES 0347−121 / BL Lacertae objects: individual: 1ES 1101−232 / BL Lacertae objects: individual: 1ES 1218+304
© A. Zech and M. Lemoine 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Blazars – active galactic nuclei (AGN) with a jetted relativistic outflow pointing towards us – are commonly classified according to their characteristic doublehumped spectral energy distribution (SED; Urry & Padovani 1995; Fossati et al. 1998). At the tip of the socalled blazar sequence (Padovani & Giommi 1995; Ghisellini et al. 2017), extreme blazars form a new population of BL Lacs (Costamante et al. 2001), designated as ‘ultrahighfrequencypeaked BL Lac objects’ (UHBLs) given the very high frequencies of their two emission peaks (Biteau et al. 2020). The subclass of ‘extremeTeV’ 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 onezone radiative models, which are generally successful in the representation of SEDs from BL Lac objects.
Although ten or so extremeTeV sources are known to date (Biteau et al. 2020), current estimates based on data from the FermiLAT γ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 extremeTeV blazars, largely simultaneous SEDs with a good multiwavelength 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 extremeTeV sources in the standard onezone 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 highfrequencypeaked 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 selfCompton (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, leptohadronic 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 ultrarelativistic 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 onezone model, adopting wellmotivated microphysical prescriptions for particle acceleration, under the assumption that electrons are coaccelerated 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 nonthermal power laws, and (2) electrons can be efficiently preheated in the transition layer of electronion 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 selfcontained 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 onezone 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 onezone 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 wellcovered 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 reacceleration 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 ep coacceleration 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}}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{e}}\phantom{\rule{0.166667em}{0ex}}\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) powerlaw population with spectral indices s_{1} (resp. s_{2}) below (resp. above) a break Lorentz factor γ_{e, b}. For such a distribution of Lorentz factors, ${\overline{\gamma}}_{\mathrm{e}}\sim \phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{e},\mathrm{min}}\phantom{\rule{0.166667em}{0ex}}{({\gamma}_{\mathrm{e},\mathrm{b}}/{\gamma}_{\mathrm{e},\mathrm{min}})}^{2{s}_{1}}$ if s_{1} < 2 < s_{2}, and ${\overline{\gamma}}_{\mathrm{e}}\sim {\gamma}_{\mathrm{e},\mathrm{min}}$ if 2 < s_{1} < s_{2}. The energy density u_{e} of the population can then be written ${u}_{\mathrm{e}}\sim {\overline{\gamma}}_{\mathrm{e}}\phantom{\rule{0.166667em}{0ex}}{n}_{\mathrm{e}}\phantom{\rule{0.166667em}{0ex}}{m}_{\mathrm{e}}{c}^{2}$.
Therefore, to convert a population of electrons (or pairs) of initial Lorentz factor γ_{0} to the above powerlaw 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_{1}.
In a (comoving) magnetic field of strength B = 10 B_{−2} mG, electrons of Lorentz factor γ_{e} = 10^{4} γ_{e 4} cool through synchrotron radiation on a length scale
$$\begin{array}{c}\hfill {l}_{\mathrm{c}}\simeq 8{\mathrm{\Gamma}}_{\mathrm{j}}\phantom{\rule{0.166667em}{0ex}}{B}_{2}^{2}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{e}\phantom{\rule{0.166667em}{0ex}}4}^{1}\phantom{\rule{0.166667em}{0ex}}\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 powerlaw 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},\mathrm{min}}^{s1}{\gamma}_{\mathrm{e},\mathrm{max}}^{3s}$ 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 β_{sh}c) with respect to the unshocked 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 RankineHugoniot conditions, or rather their relativistic generalisation (e.g., Kirk & Duffy 1999, and references therein). For a strong, subrelativistic pair shock, the postshock temperature reads ${k}_{\text{B}}{T}_{\text{e}}/{m}_{\text{e}}{c}^{2}\simeq 0.19{u}_{\text{sh}}^{2}$, while for a strong, relativistic pair shock, k_{B}T_{e}/m_{e}c^{2} ≃ 0.24u_{sh}. For reference, k_{B}T_{e}/m_{e}c^{2} ≃ 0.7 for u_{sh} = 3 in the mildly relativistic regime, close to the fully relativistic limit.
This situation changes rather dramatically in an electronion 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_{p}m_{p}c^{2} (n_{p} is understood as a proper density here), corresponding to an energyperparticle ⟨E_{p}⟩=(γ_{sh}−1)m_{p}c^{2}, 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 selfgenerated in the shock precursor (see e.g., Lemoine et al. 2019c; Vanthieghem et al., in prep.). Phenomenological models of gammaray 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 particleincell (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 Weibellike microturbulence. For reference, we define σ as the ratio of the magnetic energy density to the particle energy density in the rest frame of the unshocked (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_{B}T_{p} ≃ 0.2m_{p}c^{2} for a mildly relativistic subluminal 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_{B}T_{e} ∼ 0.1γ_{sh}m_{p}c^{2}, or a mean Lorentz factor ${\overline{\gamma}}_{\mathrm{e}}\sim 600{\gamma}_{\mathrm{sh}}$.
In the subrelativistic limit, and more generally in cases in which electron preheating is not efficient, the electrons suffer from a wellknown 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 electronion 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^{−3} (see e.g., Vanthieghem et al. 2020 and references therein). Strictly speaking, this result applies to superluminal shocks, but highly relativistic shocks are generically superluminal 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 subluminal 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 microturbulence that they themselves excite through beamplasma 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\phantom{\rule{0.166667em}{0ex}}{(4\pi {n}_{\mathrm{p}}{e}^{2}/{m}_{\mathrm{p}})}^{1/2}\sim {10}^{7}\phantom{\rule{0.166667em}{0ex}}{n}_{p\phantom{\rule{0.166667em}{0ex}}0}^{1/2}\phantom{\rule{0.166667em}{0ex}}$cm (n_{p 0} = n_{p}/1 cm^{−3}). 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^{−2}, but at least as large as the initial σ value.
The modelling of gammaray burst afterglows offers an interesting perspective on the magnitude of σ_{rad}. In those sources, the external magnetisation is significantly lower (σ ∼ 10^{−9}) than expected here, yet σ_{rad} is inferred to be as large as ∼10^{−5} − 10^{−4} (Lemoine et al. 2013; Santana et al. 2014). That effective magnetisation is generally understood as a leftover of the shockgenerated 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 RayleighTaylor 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, gammaray burst modelling provides empirical support to our model.
At moderate magnetisation, 10^{−5} ≲ σ ≲ 10^{−2}, a different currentdriven 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 extremeTeV blazars, corresponding to mildly relativistic shocks. The minimum Lorentz factor is ${\gamma}_{\mathrm{e},\mathrm{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}β_{j}c), 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 superluminal 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 ∼ Θ_{j}d 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^{−7}(n_{e}/1 cm^{−3})^{−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 superluminal and it would not arise in subluminal 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},\mathrm{max}}\lesssim \frac{{\gamma}_{\mathrm{e},\mathrm{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})^{2} and r_{L, 0} ∝ (γ_{e}/γ_{sh})/⟨B⟩. The maximum Lorentz factor arising from this limit is multiplied by m_{p}/m_{e} in an electronproton 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^{−2}, σ = 10^{−5}, γ_{sh} = 3, n_{e} = 1 cm^{−3}. Then, the scattering constraint gives γ_{e, max} ∼ 2 × 10^{5}, while the radiative and age constraints both lead to γ_{e, max} ∼ 10^{7}. 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 superluminal shocks and not subluminal 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 extremeTeV 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 subrelativistic 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; BoschRamon et al. 2012), which would instead imply a weak forward shock.
Alternatively, if an ‘obstacle’ overtakes the jet, as in the ‘blobinjet’ 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 socalled shock normal rest frame, in which the (unshocked) 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 (preshock) jet with the (oblique) shock surface, written Γ_{j<}, through Eq. (B.2): γ_{sh} = Γ_{j<}/Γ_{ns}, where Γ_{ns} 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 reconfinement 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 postshock region quite substantially relative to that in the preshock 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 δ^{4} (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 radiogalaxies, 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 preshock and postshock flows. Average is understood here over the orientation of the particle momentum in the preshock 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 onedimensional dilatation along the flow axis, $\mathrm{\Delta}lnp=\frac{1}{3}\mathrm{\Delta}ln{u}_{\mathrm{j}}$, with u_{j} the flow fourvelocity. Consider for instance a velocity profile such that Γ_{j<} decreases abruptly to Γ_{j>} through a recollimation shock, then steadily reincreases 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 subrelativistic 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 reaccelerated 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}_{>},\phantom{\rule{0.166667em}{0ex}}{\gamma}_{<})=\frac{s1}{g\phantom{\rule{0.166667em}{0ex}}{\gamma}_{<}}{\left(\frac{{\gamma}_{>}}{g\phantom{\rule{0.166667em}{0ex}}{\gamma}_{<}}\right)}^{s}\mathrm{\Theta}({\gamma}_{>}g\phantom{\rule{0.166667em}{0ex}}{\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 postshock particle distribution can be obtained through the convolution of the preshock 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}_{>}}\phantom{\rule{0.166667em}{0ex}}=\phantom{\rule{0.166667em}{0ex}}\frac{{(s1)}^{n+1}}{n!\phantom{\rule{0.166667em}{0ex}}{g}^{n}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{min}}}{\left(\frac{{\gamma}_{>}}{{g}^{n}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{min}}}\right)}^{s}ln{\left(\frac{{\gamma}_{>}}{{g}^{n}\phantom{\rule{0.166667em}{0ex}}{\gamma}_{\mathrm{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^{n}γ_{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}\mathrm{ln}\gamma $, which provides a reasonable estimate for ${\overline{\gamma}}_{\mathrm{e}}$ is found at γ = g^{n} γ_{min}exp[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^{n} γ_{e, min}, with γ_{e, min} ≃ 600γ_{sh} and the spectrum is given by Eq. (6) above.
2.5. Onezone model with e − p coacceleration in relativistic shocks
The standard onezone 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_{1} and s_{2}, and its normalisation K_{e} at γ = 1^{1}. Even with a wellsampled multiwavelength 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}, ν_{syn}L_{syn}(ν_{syn}), ν_{IC}L_{IC}(ν_{IC}), and the indices in the differential photon flux F_{ν} before and after the synchrotron peak, α_{1} and α_{2}. 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 powerlaw 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}\phantom{\rule{0.166667em}{0ex}}{e}^{\gamma /{\gamma}_{\mathrm{e},\mathrm{max}}}.\end{array}$$(7)
This simple parameterisation turns out to be well suited to characterise the SEDs of extremeTeV 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 KleinNishina 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},\mathrm{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 α_{1} ≃ 0.6 as α = (s − 1)/2 because the particles do not effectively cool through synchrotron radiation. The value of α_{2} 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\phantom{\rule{0.166667em}{0ex}}\mathrm{mG}\phantom{\rule{0.166667em}{0ex}}\frac{1+z}{\delta}\frac{{\nu}_{\mathrm{syn}}}{{10}^{17}\phantom{\rule{0.166667em}{0ex}}\mathrm{Hz}}{\left(\frac{{\gamma}_{\mathrm{e},\mathrm{max}}}{{10}^{6}}\right)}^{2}.\end{array}$$(10)
Since the above formulas were derived for a brokenpowerlaw 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 lowfrequency 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})^{2}.
The size of the homogeneous emission region is derived from the peak luminosities, following again Tavecchio et al. (1998) in the KleinNishina regime:
$$\begin{array}{c}\hfill R={\left\{2\phantom{\rule{0.166667em}{0ex}}\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}\phantom{\rule{0.166667em}{0ex}}{\left(\frac{3\delta \phantom{\rule{3.33333pt}{0ex}}{m}_{\mathrm{e}}\phantom{\rule{3.33333pt}{0ex}}{c}^{2}}{4{\gamma}_{e,max}\phantom{\rule{3.33333pt}{0ex}}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 powerlaw 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\phantom{\rule{0.166667em}{0ex}}\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 onezone leptohadronic model, this would lead to four additional free parameters for a powerlaw 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},\mathrm{min}}}^{{\gamma}_{\mathrm{p},\mathrm{max}}}\mathrm{d}\gamma \phantom{\rule{0.166667em}{0ex}}{K}_{\mathrm{p}}{\gamma}^{s+1}\phantom{\rule{0.166667em}{0ex}}{e}^{\gamma /{\gamma}_{\mathrm{p},\mathrm{max}}}\phantom{\rule{0.166667em}{0ex}}{m}_{\mathrm{p}}{c}^{2}\phantom{\rule{0.166667em}{0ex}}=\phantom{\rule{0.166667em}{0ex}}\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 preshock magnetisation. For a pure electronproton 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},\mathrm{min}}}^{{\gamma}_{\mathrm{p},\mathrm{max}}}\mathrm{d}\gamma \phantom{\rule{0.166667em}{0ex}}{\gamma}^{s}\phantom{\rule{0.166667em}{0ex}}{e}^{\gamma /{\gamma}_{\mathrm{p},\mathrm{max}}}}{{\int}_{{\gamma}_{\mathrm{e},\mathrm{min}}}^{{\gamma}_{\mathrm{e},\mathrm{max}}}\mathrm{d}\gamma \phantom{\rule{0.166667em}{0ex}}{\gamma}^{s}\phantom{\rule{0.166667em}{0ex}}{e}^{\gamma /{\gamma}_{\mathrm{e},\mathrm{max}}}}.\end{array}$$(15)
The normalisation K_{e} does not a priori fit the absolute level of the observed flux defined by ν_{syn}L_{νsyn} and ν_{IC}L_{ν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 unbeamed 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},\mathrm{max}}={\gamma}_{\mathrm{p},\mathrm{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 reaccelerated on consecutive shocks, an analogous procedure can be followed by replacing the powerlaw particle distributions with the distribution given by Eq. (6).
To enable a comparison of the energy requirements of the coacceleration model to that of standard leptonic models, the jet power is evaluated using the usual approximation (for a twosided 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 orderofmagnitude estimate.
3. Application to extremeTeV blazars
The coacceleration scenario was applied to five of the six multiwavelength SEDs from extremeTeV 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 Xray band, from NuSTAR in the hard Xray 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, extremeTeV 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^{7} s and 1ES 1218+304 at a timescale of ∼10^{5} 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 FermiLAT and H.E.S.S. spectra. Such a concave shape, if confirmed, may indicate the presence of an additional spectral component, beyond the onezone 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 nonsimultaneous 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 timedependent 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 coacceleration 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 preenergisation in some dissipation process located well upstream of the shock where electron acceleration takes place. In the latter twozone model, protons do not play any role at the shock; it would therefore apply equally well for a pure or highmultiplicity electronpositron 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; BoschRamon 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 FermiLAT 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^{4}, 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 FermiLAT 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 deabsorbed 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 coacceleration 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^{16} cm, a standard value in SSC models for highfrequencypeaked 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^{−2} is of course required to guarantee efficient acceleration.
Parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.
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: Coacceleration 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 nonrelativistic 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 selfgeneration 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 preshock 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 superluminal shocks; it would be relaxed in subluminal shocks and, in such cases, the actual preshock magnetisation could be as large as σ_{rad}. Consequently, the true physical magnetisation to be regarded is σ_{rad}, of the order of 10^{−4}–10^{−3} 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 FermiLAT flux is overpredicted 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^{−5} 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 FermiLAT 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 lowenergy 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^{3} 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 verylongbaseline 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 blobjet 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: Reacceleration 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 nonvariable 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 nonthermal emission from radio galaxies at large (kpc) distances from the core, habitually observed in the radio, optical and Xray 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 reacceleration on a second shock of a particle distribution that had already been accelerated on an initial shock. As long as cooling remains negligible, reacceleration leads to a hardening particle spectrum with a lowenergy 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 reaccelerations with a hardening power law with increasing values of ${\gamma}_{\mathrm{e},\mathrm{min}}^{(n)}\simeq {\mathit{g}}^{n}{\gamma}_{\mathrm{e},\mathrm{min}}$. If the maximum energy is limited by the scattering constraint, it increases as well by g^{n} (i.e. ${\gamma}_{\mathrm{e},\mathrm{max}}^{(n)}\simeq {\mathit{g}}^{n}{\gamma}_{\mathrm{e},\mathrm{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},\mathrm{min}}^{(n)}$ and assuming g = 2 and a power law with index s^{(0)} = 2.2 for the initial acceleration, yield indices s^{(1)} ≃ 2.1 and s^{(2)} ≃ 2.0 for the first and second reacceleration. 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 powerlaw electron distribution after n reaccelerations in a turbulent magnetic field is then given following Rybicki & Lightman (2004):
$$\begin{array}{cc}& {P}_{\mathrm{syn}}^{(n)}(\nu )\propto \phantom{\rule{0.166667em}{0ex}}{K}_{\mathrm{e}}^{(n)}\phantom{\rule{0.166667em}{0ex}}\frac{1}{{s}^{(n)}+1}\phantom{\rule{0.166667em}{0ex}}\mathrm{\Gamma}(\frac{{s}^{(n)}}{4}+\frac{19}{12})\phantom{\rule{0.166667em}{0ex}}\mathrm{\Gamma}(\frac{{s}^{(n)}}{4}\frac{1}{12})\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\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 powerlaw 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)}\phantom{\rule{0.166667em}{0ex}}[\mathrm{Hz}]\simeq 3.7\times {10}^{6}\phantom{\rule{0.166667em}{0ex}}B\phantom{\rule{0.166667em}{0ex}}[\mathrm{G}]\phantom{\rule{0.166667em}{0ex}}{{\gamma}_{\mathrm{e},\mathrm{max}}^{(n)}}^{2}$, the increase in the emitted peak power for the first reacceleration 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)})}\phantom{\rule{3.33333pt}{0ex}}\simeq \phantom{\rule{3.33333pt}{0ex}}7.1\phantom{\rule{0.166667em}{0ex}}{\left(\frac{{\gamma}_{\mathrm{e},\mathrm{max}}}{{10}^{6}}\right)}^{0.1}\phantom{\rule{3.33333pt}{0ex}}\phantom{\rule{3.33333pt}{0ex}}.\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^{3} 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 reacceleration. Following Tavecchio et al. (1998), it can be written as
$$\begin{array}{cc}& {\kappa}^{(n)}\phantom{\rule{0.166667em}{0ex}}\equiv \phantom{\rule{0.166667em}{0ex}}{\left(\frac{{L}_{\mathrm{IC}}}{{L}_{\mathrm{syn}}}\right)}^{(n)}\phantom{\rule{0.166667em}{0ex}}\propto \phantom{\rule{0.166667em}{0ex}}f({\alpha}_{1}^{(n)},{\alpha}_{2})\phantom{\rule{0.166667em}{0ex}}{\nu}_{\mathrm{s}}^{(n)}\phantom{\rule{0.166667em}{0ex}}{P}_{\mathrm{syn}}^{(n)}({\nu}_{\mathrm{s}}^{(n)})\phantom{\rule{0.166667em}{0ex}}\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\times {\left(\frac{3\phantom{\rule{0.166667em}{0ex}}{m}_{\mathrm{e}}{c}^{2}}{4\times 3.7\times {10}^{6}h}\right)}^{(3{s}^{(n)})/2}\phantom{\rule{0.166667em}{0ex}}{\left(\frac{\delta}{B\phantom{\rule{0.166667em}{0ex}}{{\gamma}_{\mathrm{e},\mathrm{max}}^{(n)}}^{3}}\right)}^{(3{s}^{(n)})/2}.\hfill \end{array}$$(19)
The ratio of the Compton dominance between the first reacceleration and the initial acceleration is then given by:
$$\begin{array}{c}\hfill \frac{{\kappa}^{(1)}}{{\kappa}^{(0)}}\phantom{\rule{3.33333pt}{0ex}}\simeq \phantom{\rule{3.33333pt}{0ex}}2.6\phantom{\rule{0.166667em}{0ex}}{\left(\frac{\delta}{50}\phantom{\rule{0.166667em}{0ex}}\frac{1\phantom{\rule{0.166667em}{0ex}}\mathrm{mG}}{B}\phantom{\rule{0.166667em}{0ex}}\frac{{10}^{6}}{{\gamma}_{\mathrm{e},\mathrm{max}}}\right)}^{0.05}.\end{array}$$(20)
The actual spectral distribution of reaccelerated 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 reacceleration 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 reacceleration 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 reaccelerations. 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 reacceleration, 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^{3} 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 onezone 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 lowenergy 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 blobinjet 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/{v}_{\perp}{\mathrm{\Gamma}}_{\text{b}}\delta $. 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 Γ_{bj} the relative blobjet Lorentz factor. On a timescale ${t}_{\text{dec}}^{\prime}$ (blob frame), the sweptup 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}}/{\mathrm{\Gamma}}_{\text{b}\text{j}}{n}_{\text{j}}{m}_{\text{p}}{c}^{2}\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_{j}m_{p}c^{2}, and Γ_{bj} 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 Poyntingflux 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 pcscale (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 superluminal 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^{−3} 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 radioloud AGN (e.g., Hada et al. 2018, and references therein), and most prominently in the most nearby radio galaxy M 87 at the ‘HST1’ radio knot (Asada & Nakamura 2012; Nakamura & Asada 2013), which appears at a deprojected 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 extremeTeV 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 abovementioned predictions from Bromberg & Levinson (2009), an orderofmagnitude 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^{2}. These estimates lead to a range of distances from the core of the order of one to several tens of parsecs.
Reacceleration 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 reacceleration.
4.2. Implications for jet power, jet content, and magnetisation
Contrary to other leptohadronic 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 leptohadronic 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 protonsynchrotron 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^{2} to 10^{3} 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^{−2}.
Studies of luminous, FSRQtype blazars show that, at least for those objects, the assumption of a pure electronproton 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 electronpositron 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 reacceleration.
The magnitude of the magnetisation σ of blazar jets is still an open question. When assuming the BlandfordZnajek process for jet launching, this implies an initially strongly Poyntingflux 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 nonstandard 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 BlandfordZnajek 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 highfrequency 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^{2} 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^{16} cm)^{3}(B_{trans}/1 G)^{4}(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 multizone 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 BlandfordPayne mechanism, a lower initial magnetisation of the jet may arise more naturally (e.g., Tzeferacos et al. 2009, for simulations in the nonrelativistic 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 frequencydependent position along the jet where synchrotron selfabsorption 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 extremeTeV 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 postshock 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 blobinjet 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 nonthermal multiwavelength emission from extremeTeV blazars to the acceleration of a protonelectron plasma on mildly relativistic shocks introduces additional constraints, arising from the microphysics in the shock region, on the standard onezone 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 coacceleration 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 extremeTeV 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 extremeTeV 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^{3} 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 reacceleration on a second shock. In the latter case, the spectral hardening and the increase in the minimum and maximum Lorentz factors arising from the reacceleration of an initial powerlaw 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
 Achterberg, A. 1990, A&A, 231, 251 [NASA ADS] [Google Scholar]
 Agarwal, A., Mohan, P., Gupta, A. C., et al. 2017, MNRAS, 469, 813 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F. A. 2000, New Astron., 5, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Araudo, A. T., BoschRamon, V., & Romero, G. E. 2010, A&A, 522, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asada, K., & Nakamura, M. 2012, ApJ, 745, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80 [Google Scholar]
 Barkov, M. V., Aharonian, F. A., & BoschRamon, V. 2010, ApJ, 724, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Begelman, M. C., & Kirk, J. G. 1990, ApJ, 353, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Biteau, J., Prandini, E., Costamante, L., et al. 2020, Nat. Astron., 4, 124 [Google Scholar]
 BoschRamon, V., Perucho, M., & Barkov, M. V. 2012, A&A, 539, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Böttcher, M., Dermer, C. D., & Finke, J. D. 2008, ApJ, 679, L9 [CrossRef] [Google Scholar]
 Bromberg, O., & Levinson, A. 2009, ApJ, 699, 1274 [CrossRef] [Google Scholar]
 Cerruti, M., Boisson, C., & Zech, A. 2013, A&A, 558, A47 [EDP Sciences] [Google Scholar]
 Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910 [Google Scholar]
 Chhotray, A., Nappo, F., Ghisellini, G., et al. 2017, MNRAS, 466, 3544 [NASA ADS] [CrossRef] [Google Scholar]
 Costamante, L. 2020, MNRAS, 491, 2771 [NASA ADS] [CrossRef] [Google Scholar]
 Costamante, L., Ghisellini, G., Giommi, P., et al. 2001, A&A, 371, 512 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Costamante, L., Bonnoli, G., Tavecchio, F., et al. 2018, MNRAS, 477, 4257 [NASA ADS] [CrossRef] [Google Scholar]
 Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, MNRAS, 485, 5105 [NASA ADS] [CrossRef] [Google Scholar]
 Dmytriiev, A., Sol, H., & Zech, A. 2021, MNRAS, 505, 2712 [NASA ADS] [CrossRef] [Google Scholar]
 Dubal, M. R., & Pantano, O. 1993, MNRAS, 261, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., & MacFadyen, A. I. 2014, ApJ, 791, L1 [CrossRef] [Google Scholar]
 Fichet de Clairfontaine, G., Meliani, Z., Zech, A., & Hervet, O. 2021, A&A, 647, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433 [Google Scholar]
 Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghisellini, G., Celotti, A., & Costamante, L. 2002, A&A, 386, 833 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011, MNRAS, 414, 2674 [NASA ADS] [CrossRef] [Google Scholar]
 Ghisellini, G., Righi, C., Costamante, L., & Tavecchio, F. 2017, MNRAS, 469, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ, 860, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Henri, G., & Pelletier, G. 1991, ApJ, 383, L7 [CrossRef] [Google Scholar]
 Hervet, O., Williams, D. A., Falcone, A. D., & Kaur, A. 2019, ApJ, 877, 26 [NASA ADS] [CrossRef] [Google Scholar]
 H.E.S.S. Collaboration (Abdalla, H., et al.) 2020, Nature, 582, 356 [Google Scholar]
 Inoue, T., Asano, K., & Ioka, K. 2011, ApJ, 734, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Katarzyński, K., Ghisellini, G., Tavecchio, F., Gracia, J., & Maraschi, L. 2006, MNRAS, 368, L52 [Google Scholar]
 Kawazura, Y., Schekochihin, A. A., Barnes, M., et al. 2020, PhRvX, 10, 041050 [NASA ADS] [Google Scholar]
 Kirk, J. G., & Duffy, P. 1999, JPhG, 25, R163 [NASA ADS] [Google Scholar]
 Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452 [NASA ADS] [Google Scholar]
 Komissarov, S. S. 2011, MmSAI, 82, 95 [NASA ADS] [Google Scholar]
 Komissarov, S. S., & Falle, S. A. E. G. 1997, MNRAS, 288, 833 [Google Scholar]
 Kumar, P., & Zhang, B. 2015, PhR, 561, 1 [NASA ADS] [Google Scholar]
 Kutkin, A. M., Sokolovsky, K. V., Lisakov, M. M., et al. 2014, MNRAS, 437, 3396 [Google Scholar]
 Lefa, E., Rieger, F. M., & Aharonian, F. 2011, ApJ, 740, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Lemoine, M. 2013, MNRAS, 428, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Lemoine, M., & Pelletier, G. 2010, MNRAS, 402, 321 [Google Scholar]
 Lemoine, M., Li, Z., & Wang, X.Y. 2013, MNRAS, 435, 3009 [NASA ADS] [CrossRef] [Google Scholar]
 Lemoine, M., Pelletier, G., Gremillet, L., & Plotnikov, I. 2014, EL, 106, 55001 [NASA ADS] [Google Scholar]
 Lemoine, M., Vanthieghem, A., Pelletier, G., & Gremillet, L. 2019a, PhRvE, 100, 033209 [NASA ADS] [Google Scholar]
 Lemoine, M., Pelletier, G., Vanthieghem, A., & Gremillet, L. 2019b, PhRvE, 100, 033210 [NASA ADS] [Google Scholar]
 Lemoine, M., Gremillet, L., Pelletier, G., & Vanthieghem, A. 2019c, PhRvL, 123, 035101 [NASA ADS] [Google Scholar]
 Martí, J. M., Perucho, M., & Gómez, J. L. 2016, ApJ, 831, 163 [Google Scholar]
 Meli, A., & Biermann, P. L. 2013, A&A, 556, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mücke, A., & Protheroe, R. J. 2001, APh, 15, 121 [Google Scholar]
 Nakamura, M., & Asada, K. 2013, ApJ, 775, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Nalewajko, K., & Sikora, M. 2009, MNRAS, 392, 1205 [CrossRef] [Google Scholar]
 O’Sullivan, S. P., & Gabuzda, D. C. 2010, ASPC, 427, 207 [Google Scholar]
 Padovani, P., & Giommi, P. 1995, ApJ, 444, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Pelletier, G., Lemoine, M., & Marcowith, A. 2009, MNRAS, 393, 587 [CrossRef] [Google Scholar]
 Pelletier, G., Gremillet, L., Vanthieghem, A., & Lemoine, M. 2019, PhRvE, 100, 013205 [NASA ADS] [Google Scholar]
 Perucho, M. 2013, EPJWC, 61, 02002 [NASA ADS] [Google Scholar]
 Petropoulou, M., & Mastichiadis, A. 2012, MNRAS, 426, 462 [NASA ADS] [CrossRef] [Google Scholar]
 Petropoulou, M., Sironi, L., Spitkovsky, A., & Giannios, D. 2019, ApJ, 880, 37 [Google Scholar]
 Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238 [Google Scholar]
 Pope, M. H., & Melrose, D. B. 1994, PASAu, 11, 175 [CrossRef] [Google Scholar]
 Potter, W. J. 2018, MNRAS, 473, 4107 [NASA ADS] [CrossRef] [Google Scholar]
 Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2009, A&A, 507, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 2004, Radiative Process in Astrophysics (Hoboken: CfA Harvard, WileyVCH) [Google Scholar]
 Saikia, P., Körding, E., & Falcke, H. 2016, MNRAS, 461, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Santana, R., Barniol, Duran R., & Kumar, P. 2014, ApJ, 785, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Sari, R., & Piran, T. 1995, ApJ, 455, L143 [Google Scholar]
 Sbarrato, T., Foschini, L., Ghisellini, G., & Tavecchio, F. 2011, AdSpR, 48, 998 [NASA ADS] [Google Scholar]
 Schneider, P. 1993, A&A, 278, 315 [NASA ADS] [Google Scholar]
 Sikora, M. 2016, Galax, 4, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Sikora, M., Madejski, G., Moderski, R., & Poutanen, J. 1997, ApJ, 484, 108 [CrossRef] [Google Scholar]
 Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103 [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54 [Google Scholar]
 Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411 [NASA ADS] [Google Scholar]
 Tavecchio, F. 2014, MNRAS, 438, 3255 [Google Scholar]
 Tavecchio, F., & Sobacchi, E. 2020, MNRAS, 491, 2198 [NASA ADS] [Google Scholar]
 Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 [Google Scholar]
 Tzeferacos, P., Ferrari, A., Mignone, A., et al. 2009, MNRAS, 400, 820 [NASA ADS] [CrossRef] [Google Scholar]
 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Vanthieghem, A., Lemoine, M., Plotnikov, I., et al. 2020, Galax, 8, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [Google Scholar]
 White, R. L. 1985, ApJ, 289, 698 [NASA ADS] [CrossRef] [Google Scholar]
 Zdziarski, A. A., Sikora, M., Pjanka, P., & Tchekhovskoy, A. 2015, MNRAS, 451, 927 [CrossRef] [Google Scholar]
 Zech, A., Cerruti, M., & Mazin, D. 2017, A& A, 602, A25 [CrossRef] [EDP Sciences] [Google Scholar]
 Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, PhRvL, 122 [Google Scholar]
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 subequipartition 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 subequipartition level by the time the electron power law is shaped.
In an electronproton 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 subrelativistic 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 nonthermal 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}}~0.10.3\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 ionelectron 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^{−3} (hence, σ_{e} ∼ 1) up to the fully relativistic regime σ_{p} ≳ 1 (σ_{e} ∼ 10^{3}). 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 nearcomplete 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 postshock flow velocity and to keep the discussion selfcontained. We do so, in particular, by deboosting to the socalled ‘shock normal frame’ ℛ_{n} (Begelman & Kirk 1990) in which the flow moves along the shock normal. This means deboosting 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 preshock from postshock quantities (indexed with _{>}). If the shock front is stationary in the source rest frame ℛ_{s}, then the shock normal frame ℛ_{n} moves at velocity β_{ns} with respect to ℛ_{s}, where
$$\begin{array}{c}\hfill {\mathit{\beta}}_{\mathbf{n}\mathbf{}\mathbf{s}}={\mathit{\beta}}_{\mathbf{j}\mathbf{<}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}({\mathit{\beta}}_{\mathbf{j}\mathbf{<}}\xb7{\mathit{n}}_{\mathbf{}\mathbf{s}}){\mathit{n}}_{\mathbf{}\mathbf{s}}\phantom{\rule{0.166667em}{0ex}}.\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 β_{ns} = β_{j<}cos(θ_{j<}−α), hence ${\mathrm{\Gamma}}_{\text{n}\text{s}}={\left[1{\beta}_{\text{j}<}^{2}{\mathrm{cos}}^{2}{\theta}_{\text{j}<}\alpha \right]}^{1/2}$.
To derive the effective velocity of the jet in this shock normal frame, we note that the component of the jet fourvelocity along the shock normal is preserved by the deboost, 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 )}\phantom{\rule{0.166667em}{0ex}}=\phantom{\rule{0.166667em}{0ex}}\frac{{\mathrm{\Gamma}}_{\mathrm{j}<}}{{\mathrm{\Gamma}}_{\mathrm{n}\mathrm{s}}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.2)
If Γ_{j<}θ_{j<}−α ≪ 1, the recollimation shock is subrelativistic in the ℛ_{n} frame, otherwise it is relativistic. The postshock 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 postshock plasma fourvelocity
$$\begin{array}{c}\hfill {\mathit{u}}_{\mathbf{j}\mathbf{>}\mathbf{}\mathbf{s}}\phantom{\rule{0.166667em}{0ex}}=\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Gamma}}_{\mathrm{j}>\mathrm{n}}{\beta}_{\mathrm{j}>\mathrm{n}}{\mathit{n}}_{\mathbf{}\mathbf{s}}\phantom{\rule{0.166667em}{0ex}}+\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Gamma}}_{\mathrm{j}>\mathrm{n}}{\mathrm{\Gamma}}_{\mathrm{n}\mathrm{s}}{\mathit{\beta}}_{\mathbf{n}\mathbf{}\mathbf{s}}\phantom{\rule{0.166667em}{0ex}},\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}}\phantom{\rule{0.166667em}{0ex}}\simeq \phantom{\rule{0.166667em}{0ex}}\frac{{\mathrm{\Gamma}}_{\mathrm{j}<}}{\sqrt{1+{\beta}_{\mathrm{j}<}^{2}{\mathrm{\Gamma}}_{\mathrm{j}<}{sin}^{2}({\theta}_{\mathrm{j}<}\alpha )}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.4)
If Γ_{j<}θ_{j<}−α ≪ 1, the shock is weak, as noted above, and the postshock plasma flows nearly along the shock surface. In the opposite limit, the recollimation shock is ultrarelativistic in ℛ_{n} (see Eq. (B.2)), the shock is strong and the postshock 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 postshock 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 postshock and preshock velocities. In the smallangle approximation, meaning θ_{j<}, θ_{j>}, α≪1, we obtain
$$\begin{array}{c}\hfill {\theta}_{\mathrm{j}>}{\theta}_{\mathrm{j}<}\phantom{\rule{0.166667em}{0ex}}\simeq \phantom{\rule{0.166667em}{0ex}}\frac{{\theta}_{\mathrm{j}<}\alpha }{\sqrt{3}}\phantom{\rule{0.166667em}{0ex}}.\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 preshock to postshock 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 comoving 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}}\phantom{\rule{0.166667em}{0ex}}{\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}}\phantom{\rule{0.166667em}{0ex}}{\delta}^{3}\phantom{\rule{0.166667em}{0ex}}{\displaystyle \int}\mathrm{d}{V}^{\prime}{j}_{{\nu}^{\prime}}^{\prime}\phantom{\rule{0.166667em}{0ex}}.\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 comoving 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}}\phantom{\rule{0.166667em}{0ex}}\frac{{\delta}^{2}}{{\mathrm{\Gamma}}_{\mathrm{b}}}\phantom{\rule{0.166667em}{0ex}}{\displaystyle \int}\mathrm{d}{V}^{\prime}{j}_{{\nu}^{\prime}}^{\prime}\phantom{\rule{0.166667em}{0ex}}\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}\phantom{\rule{0.166667em}{0ex}}=\phantom{\rule{0.166667em}{0ex}}\frac{1+z}{{d}_{\mathrm{L}}^{2}}\phantom{\rule{0.166667em}{0ex}}{\delta}^{3}\phantom{\rule{0.166667em}{0ex}}{N}_{\mathrm{b},\phantom{\rule{0.166667em}{0ex}}\mathrm{eff}}\phantom{\rule{0.166667em}{0ex}}{\displaystyle {\int}_{{V}_{1}^{\prime}}}\mathrm{d}{V}^{\prime}{j}_{{\nu}^{\prime}}^{\prime}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(C.3)
where ${V}_{1}^{\prime}$ denotes the comoving volume of one blob, of comoving 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 comoving 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}}\phantom{\rule{0.166667em}{0ex}},\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^{−1}, 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 0347121 and 1ES 1101232 for the generic model (scenario I) are shown in Figs. D.1 to D.3. For the same sources, solutions of the recollimationshock 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 reaccelerated on a second shock (scenario III) for the sources 1ES 0347121 and 1ES 1101232, which feature very hard spectra in the FermiLAT 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.
Observables and proton parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.
Fig. D.2. SSC model for the SED of 1ES 0347121 in scenario I. 
Fig. D.3. SSC model for the SED of 1ES 1101232 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 0347121 in scenario II. 
Fig. D.6. SSC model for the SED of 1ES 1101232 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 0347121 in scenario III. 
All Tables
Parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.
Observables and proton parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.
All Figures
Fig. 1. SSC model for the SED of 1ES 0229+200 in scenario I. The data points at VHE are deabsorbed 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 reacceleration 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 reaccelerations. 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 0347121 in scenario I. 

In the text 
Fig. D.3. SSC model for the SED of 1ES 1101232 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 0347121 in scenario II. 

In the text 
Fig. D.6. SSC model for the SED of 1ES 1101232 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 0347121 in scenario III. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.