Electron-proton co-acceleration on relativistic shocks in extreme-TeV blazars

The multi-wavelength emission from a newly identified population of `extreme-TeV' blazars, with Compton peak frequencies around 1 TeV, is difficult to interpret with standard one-zone emission models. Large values of the minimum electron Lorentz factor and quite low magnetisation values seem to be required. We propose a scenario where protons and electrons are co-accelerated 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. The shock co-acceleration scenario provides additional constraints on the set of parameters of a standard one-zone lepto-hadronic 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 extreme-TeV sources we have examined, re-acceleration on a second shock appears needed for those objects with the hardest gamma-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 self-consistent representations were found for the most prominent representatives of this new blazar class.


Introduction
Blazars -active galactic nuclei (AGN) with a jetted relativistic outflow pointing towards us -are commonly classified according to their characteristic double-humped spectral energy distribution (SED; Urry & Padovani 1995;Fossati et al. 1998). At the tip of the so-called blazar sequence (Padovani & Giommi 1995;Ghisellini et al. 2017), extreme blazars form a new population of BL Lacs (Costamante et al. 2001), designated as 'ultra-highfrequency-peaked BL Lac objects' (UHBLs) given the very high frequencies of their two emission peaks (Biteau et al. 2020). The sub-class of 'extreme-TeV' sources presents rather unusual properties, with the higher energy peak beyond ∼ 1 TeV and a hard spectrum in the γ-ray range. These objects still defy a theoretical interpretation in terms of the commonly used one-zone radiative models, which are generally successful in the representation of SEDs from BL Lac objects.
Although ten or so extreme-TeV sources are known to date (Biteau et al. 2020), current estimates based on data from the Fermi-LAT γ-ray telescope predict about five times more candidate sources that need to be confirmed with future instruments in the TeV range (Costamante 2020). For five of the known extreme-TeV blazars, largely simultaneous SEDs with a good multi-wavelength coverage from the infrared to the TeV band have become available, and they have already provided important clues on the physics of dissipation in these sources (Costamante et al. 2018).
Arriving at a satisfactory description of extreme-TeV sources in the standard one-zone leptonic framework seemingly requires two essential ingredients (e.g. Tavecchio, Maraschi, & Ghisellini 1998;Ghisellini, Celotti, & Costamante 2002;Katarzyński et al. 2006;Costamante et al. 2018): (1) lower magnetic fields ( 10 mG) than for more common high-frequency-peaked BL Lac objects (HBLs), corresponding to magnetic field energy densities well below equipartition, and (2) a peculiar lepton distribution with most of the energy carried by particles of a large Lorentz factor, γ e ∼ 10 3 − 10 4 . The former requirement avoids an excessive softening of the γ-ray spectrum by synchrotron cooling and is necessary to reproduce the large separation between the synchrotron and synchrotron self-Compton (SSC) peak frequencies, while the latter is needed to accommodate the steep γ-ray spectrum and high peak frequency. Alternative models have been devised (e.g. Aharonian 2000;Böttcher, Dermer, & Finke 2008;Lefa, Rieger, & Aharonian 2011;Tavecchio 2014;Cerruti et al. 2015;Chhotray et al. 2017;Tavecchio & Sobacchi 2020), but they all come at the price of greater complexity and most often invoke other extreme parameter values, which often lack experimental or physical justification. For instance, lepto-hadronic models can provide reasonable SEDs with less extreme values of γ e , but at the expense of a high jet power, mostly carried by the accompanying population of ultra-relativistic protons or the magnetic field (e.g. Cerruti et al. 2015).
Article number, page 1 of 18 A&A proofs: manuscript no. zech_lemoine_2021_final In the present study we take at face value the low magnetisation inferred from the SEDs of extreme blazars and propose revisiting the one-zone model, adopting well-motivated microphysical prescriptions for particle acceleration, under the assumption that electrons are co-accelerated with protons in relativistic internal or recollimation shocks. Two key motivations of our study are the following: (1) at low magnetisation, relativistic shock acceleration is expected to dominate the physics of dissipation and of the generation of non-thermal power laws, and (2) electrons can be efficiently preheated in the transition layer of electron-ion relativistic shocks, up to a fraction of the thermal energy of shocked ions, implying large minimum Lorentz factors for the electron population, as seemingly required by phenomenology (e.g. Vanthieghem et al. 2020, and references therein). We further aim at making our approach as self-contained as possible, by fixing the microphysical parameters describing the acceleration process from first principle studies of relativistic shock models. Our main objective is to provide a critical discussion of such models to reproduce the SEDs of extreme blazars, in a modern context. This contrasts with generic one-zone models, which instead parameterise the accelerated particle population and adjust the parameters to the observed SEDs.
We organise the discussion as follows. In Sect. 2 we discuss existing constraints on the microphysical parameters of shock acceleration. The standard one-zone SSC model is then extended to include a population of protons, and physical constraints from the considered acceleration mechanism are translated into constraints on the parameters that describe the two particle populations. Section 3 then applies the resulting model to the set of well-covered SEDs. In a first scenario, generic solutions with a minimum number of free parameters are presented. A second scenario applies realistic parameters for acceleration on a mildly relativistic internal or recollimation shock, while a third scenario accounts for the possibility of re-acceleration at multiple shock waves, which appears necessary to reproduce the hard γray spectra of certain objects. A critical discussion of this model, as well as the direct consequences for our understanding of the jets of extreme blazars, is given in Sect. 4.

A large γ e
We define the Lorentz factor γ e as the average over the population, with n e = dγ e dn e /dγ e the electron number density. Blazar SEDs are generically well represented by radiative leptonic emission from a (broken) power-law population with spectral indices s 1 (resp. s 2 ) below (resp. above) a break Lorentz factor γ e,b . For such a distribution of Lorentz factors, γ e ∼ γ e,min γ e,b /γ e,min 2−s 1 if s 1 < 2 < s 2 , and γ e ∼ γ e,min if 2 < s 1 < s 2 . The energy density u e of the population can then be written u e ∼ γ e n e m e c 2 .
Therefore, to convert a population of electrons (or pairs) of initial Lorentz factor γ 0 to the above power-law spectrum necessitates the dissipation of an energy reservoir that is γ e /γ 0 times larger than the initial particle energy. In this sense, the large inferred value of γ 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), γ e ∼ 10 3 − 10 4 arises as a generic prediction independently of the value of s 1 .
In a (co-moving) 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 as measured in the source rest frame, and expressed in terms of Γ j the jet Lorentz factor. Hence, one cannot a priori exclude that electrons gained their large average Lorentz factors in some primary dissipation region, located well upstream of the acceleration region where the power-law population is generated. Such a scenario is nevertheless constrained by the apparent lack of radiation associated with that putative primary dissipation process, as well as by the shortened cooling length in the primary dissipation region with enhanced magnetic strength. In the absence of a definite model for the evolution of the magnetic field strength, electron distribution etc. along the jet, it is difficult to make this statement more quantitative as the synchrotron luminosity scales with the combination Γ 4 j n e R 3 u B γ s−1 e,min γ 3−s e,max 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.

Electron heating in relativistic shocks
Consider a shock front, moving at Lorentz factor γ sh (and velocity β sh c) with respect to the un-shocked plasma. In the context of extreme blazars, we are mostly interested in the case of mildly relativistic or truly relativistic shocks (i.e. u sh ≡ γ sh β sh 1). We also assume here a shock of normal incidence, meaning that the upstream flow enters the shock front along the shock normal, whose direction is perpendicular to the shock front. We generalise this to the case of oblique shocks in Sect. 2.4.2.
If the composition is purely leptonic, the electrons are energised through shock crossing according to the standard Rankine-Hugoniot conditions, or rather their relativistic generalisation (e.g. Kirk & Duffy 1999, and references therein). For a strong, sub-relativistic pair shock, the post-shock temperature reads k B T e /m e c 2 ≃ 0.19u 2 sh , 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 electron-ion shock, assuming for simplicity an equal number of ions and electrons. Most of the energy that enters the shock front is now carried by the ions, but a fraction of it is given to the electrons, which are thereby preheated up to a fraction of equipartition. More specifically, in the reference frame in which the shock front lies at rest, the kinetic energy density in (cold) protons is e p = γ sh (γ sh − 1) n p m p c 2 (n p is understood as a proper density here), corresponding to an energy-per-particle 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 self-generated in the shock precursor (see e.g. Lemoine et al. 2019a;Vanthieghem et al. 2021). Phenomenological models of gamma-ray burst afterglows provide a nice empirical confirmation of this energy transfer between ions and electrons, as the electron energy fraction is almost always found to be within a factor of a few from the ion energy fraction (e.g. Kumar & Zhang 2015, and references therein).
In the absence of a definitive model, we can rely on particle-in-cell (PIC) simulations to provide estimates for the average electron Lorentz factor. Generally speaking, the preheating appears more efficient at lower magnetisations σ (Sironi, Spitkovsky, & Arons 2013), likely because of a larger precursor and because of the presence of intense Weibel-like microturbulence. For reference, we define σ as the ratio of the magnetic energy density to the particle energy density in the rest frame of the un-shocked (ambient) plasma. Simulations also report a more efficient preheating in the fully relativistic regime -meaning γ sh 3 − 5 -than in the mildly relativistic regime γ sh ∼ 1 − 3. For example, Crumley et al. (2019) measure T e ≃ 0.3T p and k B T p ≃ 0.2m p c 2 for a mildly relativistic sub-luminal shock of velocity β sh = 0.8, corresponding to γ sh = 1.7, at magnetisation σ ∼ 0.007 (here defined downstream). This electronic temperature corresponds to a mean Lorentz factor γ e ∼ 300. At larger shock velocities, and smaller magnetisations, Sironi & Spitkovsky (2011) and Sironi, Spitkovsky, & Arons (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 γ e ∼ 600γ sh .
In the sub-relativistic limit, and more generally in cases in which electron preheating is not efficient, the electrons suffer from a well-known injection problem into the Fermi acceleration cycle. In such cases, the energy density of the suprathermal tail of accelerated electrons represents a modest fraction ǫ e of the energy that is incoming into the shock, for example ǫ e ≃ 0.001 for the case with β sh = 0.8 discussed in Crumley et al. (2019). Correspondingly, the minimum Lorentz factor of the suprathermal electron power law -that is, the Lorentz factor at which the power law emerges out of the thermal Maxwellian -is substantially larger than γ 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 γ 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 γ e , but also a hard spectrum with index s ∼ 2. In the following, we thus adopt ǫ e = 0.1 and γ e ∼ 600γ sh as fiducial values, corresponding to relativistic, weakly magnetised electron-ion shocks. We also assume one electron per ion, as larger electron multiplicities are expected to lead to lower mean energy per particle, although this case has not been explicitly probed by numerical experiments.

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 super-luminal shocks, but highly relativistic shocks are generically super-luminal as a consequence of the Lorentz boost (a factor γ sh ) of the perpendicular magnetic field components in the shock rest frame. In the mildly relativistic limit, however, such effects disappear and sub-luminal configurations become about as likely. Consequently, particle acceleration can become efficient, even at relatively high magnetisations (see e.g. Crumley et al. 2019) In the regime of low magnetisation that we are interested in, the size of the precursor is governed by the scattering of particles in the electromagnetic micro-turbulence that they themselves excite through beam-plasma instabilities (e.g. Plotnikov, Grassi, & Grech 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;Pelletier et al. 2019;Lemoine et al. 2019b,c). 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 λ p ≡ c 4πn p e 2 /m p −1/2 ∼ 10 7 n 1/2 p 0 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 gamma-ray 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, Li, & Wang 2013;Santana, Barniol Duran, & Kumar 2014). That effective magnetisation is generally understood as a left-over 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, Asano, & Ioka 2011) or a Rayleigh-Taylor instability at the contact discontinuity (Duffell & MacFadyen 2014). The value of σ rad that we use in our models is comparable to that above; in that sense, gammaray burst modelling provides empirical support to our model. At moderate magnetisation, 10 −5 σ 10 −2 , a different current-driven instability is expected to contribute to the generation of the magnetised microturbulence (Lemoine et al. 2014). The consequences of this instability, which is driven by the perpendicular current carried by the accelerated particles around the mean background magnetic field, have not been examined to the extent needed. In a first approximation, we can consider that the turbulence will have the same general characteristics as that generated by the current filamentation instability, that is, a strength ǫ B ∼ 0.01 and a length scale ∼ O(λ 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, Article number, page 3 of 18 A&A proofs: manuscript no. zech_lemoine_2021_final the spectral index is expected to take values s ≃ 2 − 2.3. For definiteness, we assume s ≃ 2.2 in our application to extreme-TeV blazars, corresponding to mildly relativistic shocks. The minimum Lorentz factor is γ e,min ≃ γ e ≃ 600γ 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 super-luminal shocks) t scatt r L,0 , where r L,0 represents the particle gyroradius in the coherent (background) magnetic field. We note that the limit associated with lateral escape through the jet boundary, viz. t acc r, is comparable to the age limit, if r ∼ Θ 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, Lemoine, & Marcowith 2009;Lemoine & Pelletier 2010). This constraint thus assumes the magnetic field to be super-luminal 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: 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 electron-proton plasma. (See Vanthieghem et al. (2020) for more details.) In the models that we adjust to the blazar SEDs further below, we find that this third constraint provides the limiting factor, although the radiative constraint and the age constraints do not lag by much. To put in some numbers, assume ǫ B = 10 −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 super-luminal shocks and not sub-luminal ones.
We also modelled the radiation from protons, although we find that it is always negligible for the model parameters that we derive. The spectrum characterising the protons can be similarly defined by a minimum Lorentz factor, γ p,min ∼ γ sh , a maximal Lorentz factor determined by the conjunction of the age and scat-tering constraints above and the same spectral index as for the electrons.

Shocks of normal incidence
The interaction of a jet with an obstacle generally triggers a double shock configuration, including a forward and a reverse shock. Depending on the relative momentum fluxes carried by the jet and the obstacle, either or both of these shocks can be strong or weak (Sari & Piran 1995). Given that particles release their radiation in the downstream rest frame of the shock where they are accelerated, the velocity of that downstream frame sets the Doppler factor that modulates the radiation. For extreme-TeV blazars, the large inferred Doppler boosting indicates that this downstream frame moves at large relativistic velocities towards the observer in the source frame. For shocks of normal incidence, this implies in turn that the shock itself moves at large velocities in the source frame, because the downstream rest frame moves at sub-relativistic velocities relative the shock.
Consequently, if the jet overtakes an obstacle, we need to assume that the jet ram pressure far exceeds that of the obstacle, so that the Doppler boosting of the emission is given by that of the jet itself. The shock Lorentz factor γ sh corresponds here to the relative Lorentz factor between the shock front and the obstacle, and it is of the order of Γ j if the obstacle moves at slow velocities in the source rest frame. This may give rise to shocks of large Lorentz factors. However, obstacles can penetrate a jet only if their ram pressure exceeds that of the jet (e.g. Araudo, Bosch-Ramon, & Romero 2010; Barkov, Aharonian, & Bosch-Ramon 2010; Bosch-Ramon, Perucho, & Barkov 2012), which would instead imply a weak forward shock.
Alternatively, if an 'obstacle' overtakes the jet, as in the 'blob-in-jet' scenario, the Doppler factor will be at least as large as that of the jet itself. If the blob ram pressure exceeds that of the jet, the Lorentz factor of the emission region is given by that of the blob, Γ b , while the shock Lorentz factor γ sh is set as before by the relative Lorentz factor between the blob and the jet. From Mildly relativistic shocks can be attained for sufficiently large ratios of Γ b to Γ j . We consider such a configuration in Sect. 3.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í, Perucho, & Gómez 2016, and references therein). The general properties of these recollimation shocks have been explored in various studies (e.g. Dubal & Pantano 1993;Komissarov & Falle 1997;Nalewajko & Sikora 2009;Bromberg & Levinson 2009). We summarise in Appendix B the characteristics of interest to us, in particular the jump conditions, and we establish here the correspondence between these properties and the above constraints on particle acceleration. In brief, these constraints remain applicable to oblique shocks once they are expressed in the so-called shock normal rest frame, in which the (un-shocked) flow impinges on the shock at normal incidence (Begelman & Kirk 1990). The Lorentz factor γ sh , which corresponds to the Lorentz factor of the shock with respect to the upstream plasma in this shock normal frame, is writ-ten Γ j<|n in Appendix B. It can be related to the relative Lorentz factor of the (pre-shock) jet with the (oblique) shock surface, written Γ j< , through Eq. (B.2): γ sh = Γ j< /Γ n|s , where Γ n|s represents the Lorentz factor associated with the boost to the shock normal frame. The main result here can be phrased in terms of the angle θ j< − α between the flow incidence and the shock surface (θ j< and α, respectively, represent the angles of the flow and the shock surface with respect to the jet axis). Given that the jet opening angle is of the order of 1/Γ j , one expects θ j< − α ∼ 1/Γ j in a simple re-confinement scenario (Nalewajko & Sikora 2009), and we therefore write θ j< − α ∼ κ/Γ j< . Then, we obtain γ sh ≃ √ 1 + κ 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. 2020).
Similarly, we can write the Lorentz factor of the flow past the shock as Γ j> ≃ Γ n|s ≃ Γ j< / √ 1 + κ 2 . This Lorentz factor is that which controls the amount of Doppler boosting from the emission region, provided the recollimation shock is stationary in the source frame. We note that the flow is refracted by an angle of the order of θ j< − α in the source frame as it crosses the shock, which may therefore affect the overall Doppler boosting from the post-shock region quite substantially relative to that in the pre-shock region. We do not consider this effect in the models that follow, for the sake of simplicity.
In our models below, we seek to describe the radiation emitted by an element of plasma crossing such a shock, or a series of shocks. In this description, we can model the plasma element as a blob containing a fixed number of particles, in which particles get accelerated at some stage then radiate. Even though we are dealing with a standing shock, which provides a steady emission pattern, the corresponding Doppler boosting of the radiation from the blob, as seen by the distant observer, is a factor δ 4 (c.f. Sikora et al. 1997), with δ the bulk Doppler factor. For the sake of completeness, we provide more details on this issue in Appendix C.

Interaction with multiple shocks
As recollimation (or, more generally, standing) shocks may come in series of shock fronts in the jets of radio-galaxies, we must consider the possibility that the radiating blob undergoes multiple episodes of shock acceleration, which may substantially modify the spectrum of accelerated particles (White 1985;Achterberg 1990;Schneider 1993;Pope & Melrose 1994;Meli & Biermann 2013). Particles are heated and accelerated in those shocks, but they also suffer adiabatic cooling (and possibly, radiative cooling) in the rarefaction regions that separate them.
As a particle crosses a shock front, its energy in the local frame of rest of the plasma is increased on average by the relative Lorentz factor between pre-shock and post-shock flows. Average is understood here over the orientation of the particle momentum in the pre-shock flow. For an oblique shock, the relative Lorentz factor Γ rel can be written in terms of normal shock frame quantities Γ rel ≃ 0.7Γ j<|n , where the prefactor holds for mildly relativistic shocks with Γ j<|n 2 − 3. In terms of source frame quantities, Γ rel ≃ 0.5Γ j< /Γ j> (in analogy with Eq. 4).
Adiabatic cooling in the rarefaction regions implies, for onedimensional dilatation along the flow axis, ∆ ln p = − 1 3 ∆ ln u j , with u j the flow four-velocity. Consider for instance a velocity profile such that Γ j< decreases abruptly to Γ j> through a rec-ollimation shock, then steadily re-increases to its initial value Γ j< through rarefaction. Overall, the momentum of the particle increases by g ∼ Γ 2/3 rel between a time defined as immediately before the crossing of the first shock and that immediately before the crossing of the second shock. The difference with the sub-relativistic regime studied in the above references should be noted : there the particle gains little energy through shock crossing, but cools efficiently in the rarefaction region, leading to an overall energy loss in between two consecutive shock crossings.
Beyond this sequence of heating and cooling, which shifts the energy spectrum in momentum space while preserving its shape, the particles are also re-accelerated at the shock. This now modifies the spectral shape. This effect can be described, in a first approximation, through the effective propagator 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, g ≃ Γ 2/3 rel as discussed above. The post-shock particle distribution can be obtained through the convolution of the pre-shock distribution with this propagator, which gives, after n shock crossings (White 1985;Achterberg 1990;Schneider 1993) 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 n−th shock. The hardening is stronger at momenta close to the effective injection momentum than at high energy (HE). In particular, the maximum of dN (n) > /d ln γ, which provides a reasonable estimate for γ e is found at γ = g n γ min exp [n/(s − 1)], which indicates a rapid rise of the effective γ 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.

One-zone model with e − p co-acceleration in relativistic shocks
The standard one-zone SSC model has usually nine free parameters: the strength of the uniform magnetic field B, the size of the spherical emission region R, its Doppler factor δ, the minimum, maximum and break Lorentz factors of the electron distribution γ e,min , γ e,max , γ e,br , the two indices of the electron distribution before and after the break s 1 and s 2 , and its normalisation K e at γ = 1 1 . Even with a well-sampled multi-wavelength dataset, this model remains highly degenerate. Tavecchio, Maraschi, & Ghisellini (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, Boisson, & Zech (2013).
The constraints on the magnetisation, on γ e,min , γ e,max and on the spectral index, which we derived above from the microphysics of relativistic shock acceleration, remove most of the inherent degeneracy of the standard model. We demonstrate this by considering a power-law electron distribution with an exponential cutoff of the form This simple parameterisation turns out to be well suited to characterise the SEDs of extreme-TeV blazars.
For an initially chosen value of δ, first the maximum electron Lorentz factor γ e,max can be determined from the IC peak frequency ν IC . In the Klein-Nishina regime, which is relevant for the datasets we want to investigate in the following section, this relation is given following Tavecchio, Maraschi, & Ghisellini (1998): with g(α 1 , α 2 ) = exp 1 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 (10) Since the above formulas were derived for a broken-powerlaw 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, Maraschi, & Ghisellini (1998) in the Klein-Nishina regime: 3δ m e c 2 4γ e,max hν syn with We may then verify our assumption of a single power-law distribution without cooling break, by comparing the frequency at which synchrotron cooling becomes relevant, namely versus γ e,max . Given the low magnetisation required for efficient shock acceleration, the assumption holds in our application. IC cooling is not significant either in the absence of strong Compton dominance.
Hence, the only missing parameter is the normalisation of the electron distribution K e . Following our scenario outlined in the previous section, a population of protons is added to the model. In a standard one-zone lepto-hadronic model, this would lead to four additional free parameters for a power-law distribution (e.g. Cerruti, Boisson, & Zech 2013;Zech, Cerruti, & Mazin 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 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 electron-proton plasma with equal particle number densities n e = n p , we then find 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 un-beamed sources (Tavecchio, Maraschi, & Ghisellini 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 γ p,max = γ p,min / √ σ, in line with the scalings adopted for the electrons, up to preheating in the shock transition. When considering a scenario where particles are re-accelerated on consecutive shocks, an analogous procedure can be followed by replacing the power-law particle distributions with the distribution given by Eq. (6).
To enable a comparison of the energy requirements of the co-acceleration model to that of standard leptonic models, the jet power is evaluated using the usual approximation (for a twosided jet): where u γ represents the energy density contained in the SSC radiation field. We assume a value for the bulk Lorentz factor Γ = 0.5 δ, corresponding to a situation where the jet is closely aligned with the line of sight. The required luminosity L j needs to be compared to the theoretically available power from the central engine, for which the Eddington luminosity of the central black hole can be used as an order-of-magnitude estimate.

Application to extreme-TeV blazars
The co-acceleration scenario was applied to five of the six multi-wavelength SEDs from extreme-TeV blazars presented by Costamante et al. (2018) for the sources 1ES 0229+200, 1ES 1101-232, 1ES 0347-121, RGB J0710+591, 1ES 1218+304, and 1ES 0414+009. These datasets were chosen due to the dense and homogeneous wavelength coverage including simultaneous data from Swift in the UV and soft X-ray band, from NuStar in the hard X-ray band and Fermi in the HE γ-ray band, together with contemporaneous data at very high energy (VHE) from H.E.S.S. and VERITAS, as well as archival data from WISE in the infrared. It should be noted that, although the VHE data are not strictly simultaneous, extreme-TeV blazars are known for their lack of variability in the VHE range. Some level of variability has so far only been observed for 1ES 0229+200 at a timescale of ∼ 10 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 Fermi-LAT and H.E.S.S. spectra. Such a concave shape, if confirmed, may indicate the presence of an additional spectral component, beyond the one-zone description we are providing here. It is however plausible, as assumed by Costamante et al. (2018), that a temporal variation in the spectral shape in the HE band is responsible for this apparent upturn between the non-simultaneous datasets. Given these uncertainties, this source has thus been excluded from our study, reducing the sample to five sources.
As a first step, the values of the model parameters were determined using the analytical equations given in Sect. 2.5. A high Doppler factor of δ = 50 was chosen initially following the results of Costamante et al. (2018) and adjusted to lower or higher values if a good representation of the data could not be achieved. The estimated parameter values were used as a starting point for modelling with a numerical code that treats the radiative transfer and all relevant leptonic and hadronic emission processes (Cerruti et al. 2015;Zech, Cerruti, & Mazin 2017). The contribution from the host galaxy was reproduced with a template from Silva et al. (1998), where the normalisation was adjusted to fit the infrared and optical data points. We have verified, analytically and with the time-dependent code by Dmytriiev, Sol, & Zech (2021), that in the resulting solutions, radiative cooling of the emitting particle population can be safely neglected, given the low magnetisation that is inherent to the shock acceleration scenario.
The co-acceleration model was declined into three scenarios to find the physically most meaningful representations of the selected SEDs. In a first scenario, we study a generic version with a large value of γ e,min , mimicking electron preheating at a shock of large Lorentz factor, or in some unspecified primary dissipation region. In a second scenario, we study the case of a lower value of γ e,min , as discussed earlier in the context of preheating at mildly relativistic internal or recollimation shocks. Finally, in a third scenario, we investigate the possibility of reacceleration at multiple shocks, focussing on the case of a single reacceleration at a secondary shock, sufficient to provide a good representation for the datasets under study, as will be seen.
In all cases, synchrotron radiation from the protons peaked at radio frequencies and was several orders of magnitude lower than the emission from electrons. It can thus be completely neglected when studying the SEDs. Only the kinetic energy of the protons needs to be considered when discussing the energetics of the source.
3.1. Scenario I: Generic shock model with large γ e,min We thus consider here a model with a large value of γ e,min . In light of our discussion in Sect. 2, such a large value can either represent the consequence of electron preheating at a shock of large Lorentz factor or some pre-energisation in some dissipation process located well upstream of the shock where electron acceleration takes place. In the latter two-zone model, protons do not play any role at the shock; it would therefore apply equally well for a pure or high-multiplicity electron-positron composition.
In the former, a large value of the shock Lorentz factor (e.g. ∼ δ), means that the jet encounters an obstacle moving at small velocity in the source rest frame (see Sect. 2.4.1). However, this requires the jet ram pressure to exceed that of the obstacle in the source frame, which is precisely opposite to the condition that allows the obstacle to penetrate the jet (e.g. Araudo, Bosch-Ramon, & Romero 2010; Barkov, Aharonian, & Bosch-Ramon 2010; Bosch-Ramon, Perucho, & Barkov 2012). Consequently, none of the above possibilities currently receives clear justification. For this reason, we regard this first scenario as a case study. To fix the parameters, we assume that γ e,min is set by the shock Lorentz factor, in accordance with Sect. 2.3, and that this shock Lorentz factor is set equal to the Doppler factor δ, then slightly adjusted by constraints from the slope of the Fermi-LAT spectrum where necessary. This turned out to be necessary for one source only, 1ES 1218+304 (see below).
This choice leads to very high values of γ e,min , of a few 10 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. Fig. D.1  to D.3) if, in addition to adjusting δ to values between 30 and 50, one also allows for an adjustment of the size of the emission region R, typically by a factor of ∼ 1.5. In the case of 1ES 1218+304, the Fermi-LAT spectrum is flatter than for the other sources, requiring a decrease in γ e,min by an order of magnitude, in addition to adjustments of δ and R (cf. Fig. 2).
Allowing for these adjustments, satisfactory solutions can be found for all five SEDs. The full application of the constraints of the shock co-acceleration model leads to unique solutions in this generic scenario. As can be seen from Table 1, the parameter values are very similar for the different sources: the Doppler factor of the emission region is large (up to δ = 60); the emission region has a size of a few 10 16 cm, a standard value in SSC models for high-frequency-peaked BL Lac objects (HBLs); the magnetic field strength, in the milligauss range, and the magnetisation are very small, as has been found before for these extreme sources. In our framework, a low magnetisation σ ≪ 10 −2 is of course required to guarantee efficient acceleration.
In certain cases, for example the model of 1ES 1218+304 shown in Fig. 2, a mismatch between the model and the SED appears in the ultraviolet and infrared ranges. This may be explained with an additional synchrotron emission component from the extended jet, which is not represented in our models.

Scenario II: Co-acceleration on mildly relativistic shocks
In this second scenario, we consider the more physically grounded case of electron preheating and acceleration at internal (Sect. 2.4.1) or recollimation shocks (Sect. 2.4.2). In both cases, far smaller values of the shock Lorentz factor in the shock normal frame Γ j<|n should be assumed than those applied in the generic scenario I.
For acceleration on an internal shock caused by a plasma blob with Lorentz factor Γ b traversing the jet, the resulting shock Lorentz factor γ sh ≡ Γ j<|n is given by Eq. 4. Assuming a Doppler factor of δ for the observed emission from the blob, a mildly relativistic shock with Γ j<|n ≃ 3 would result for a jet Lorentz factor Γ j ≃ δ/12, when assuming a closely aligned jet, corresponding to Γ b ≃ δ/2, such that Γ j would be of the order of a few. Such a combination of a highly relativistic emission component responsible for HE emission and a mildly relativistic jet is frequently proposed to explain the discrepancy in observed bulk Lorentz factors in the radio and HE bands (Ghisellini, Tavecchio, & Chiaberge 2005;Henri & Pelletier 1991;Sol, Pelletier, & Asseo 1989). In this case, shock Lorentz factors much higher than a few can only be achieved for unrealistically high blob Lorentz factors or for non-relativistic jets.
Considering the second case of acceleration on recollimation shocks, since the bulk Lorentz factor of the downstream, radiating plasma is Γ j>|s ≥ δ/2, the shock Lorentz factor must verify Γ j<|n ≃ Γ j<|s /Γ j>|s 2Γ j<|s /δ (cf. Appendix B). If the bulk Lorentz factor of the plasma flow upstream of the recollimation shock Γ j<|s is of a similar strength as the typical bulk Doppler factors derived from SED modelling, the shock Lorentz factor Γ j<|n should be close to unity. The default value of Γ j<|n = 3 that we impose in this second scenario for both cases to better accommodate the datasets is already relatively high, but it should still be acceptable in this case.
With this assumption, the minimum electron Lorentz factor is fixed to γ e,min ≃ 600 Γ j<|n ≃ 1800. The corresponding model fits imply different values of the physical parameters, following the procedure outlined in Sect. 2.5. Retaining bulk Doppler factors around δ ≈ 50, as for the solutions in scenario I, the decrease in γ e,min for a fixed γ e,max leads to a quadratic decrease in the inferred upstream magnetisation σ (cf. Eq. (3)). Since the particle density is inversely proportional to the magnetisation downstream of the shock [cf. Eq. (14)], the condition σ rad = σ needs to be abandoned to avoid a strong increase in the particle density, which would lead to a Compton dominance in the SED that is not observed in extreme blazars.
An increase in the magnetisation in the region downstream of the shock by one or two orders of magnitude is not unexpected, due to the self-generation of magnetic turbulence at the shock. As discussed earlier, shock microphysics only imposes the condition σ ≤ σ rad 0.01, with 0.01 representing the typical effective magnetisation at the shock, σ rad the magnetisation in the radiation region and σ the pre-shock magnetisation. We also stress that the values of σ that are provided in Tab. 1 are only indicative, in the sense that they are derived from the respective values of γ e,min and γ e,max under the assumption that the maximal energy is set by the scattering constraint. As discussed earlier, this constraint applies formally to super-luminal shocks; it would be relaxed in sub-luminal shocks and, in such cases, the actual pre-shock magnetisation could be as large as σ rad . Consequently, the true physical magnetisation to be regarded is σ rad , of the order of 10 −4 to 10 −3 for the various sources in this model.
As can be seen from Figures 3, D.4 and D.5, satisfactory representations of the SEDs can be found within this shock scenario for 1ES 1218+304, RGB J0710+591 and 1ES 0347-121. While for the first two, the representation of the SED is of comparable quality to the one in scenario I, for 1ES 0347-121 the Fermi-LAT flux is over-predicted for the lowest flux point, but given the uncertainties in the Fermi points, this solution can still be regarded as acceptable. It should be noted here that, as far as the synchrotron and SSC peak flux levels are concerned, the residual freedom in σ rad introduces a degeneracy in the sense that similar solutions can be found by increasing δ and decreasing σ rad or vice versa.
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 Fermi-LAT points, the infrared flux is overproduced by the model. If γ e,min remains fixed at its default value, the only way around this problem is to increase δ to a value well above 100, in order to shift the 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 Fig. 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 very-long-baseline interferometry (VLBI) and on models of relativistic jets, are of the order of a few to a few tens (Pushkarev et al. 2009;Saikia, Körding, & Falcke 2016).
For acceleration on internal shocks caused by blob-jet interaction, the required values of Γ j<|n can be accommodated with jet Lorentz factors of a few, if the jet is not too closely aligned with the line of sight (i.e. for δ Γ b ).

Scenario III: Re-acceleration on a second shock
Despite the difficulties the recollimation shock scenario has in accounting for the SEDs with the hardest γ-ray spectra, it remains attractive as a natural explanation for continuous 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 non-thermal emission from radio galaxies at large (kpc) distances from the core, habitually observed in the radio, optical and X-ray bands, and recently even in VHE γrays (H. E. S. S. Collaboration et al. 2020). For simplicity, and to keep the model as economical as possible, we consider the simplest generalisation of scenario II, namely re-acceleration on a second shock of a particle distribution that had already been accelerated on an initial shock. As long as cooling remains negligible, re-acceleration leads to a hardening particle spectrum with a low-energy turnover at increasingly higher energies above the initial value of γ e,min , as detailed in Sect. 2.4.3.
The effect on the SED can be roughly estimated by approximating the electron distribution after n re-accelerations with a hardening power law with increasing values of γ (n) e,min ≃ g n γ e,min . If the maximum energy is limited by the scattering constraint, it increases as well by g n (i.e. γ (n) e,max ≃ g n γ e,max ), provided the magnetisation remains comparable from shock to shock. It cannot, of course, exceed the synchrotron or age limit. Fits of a power law to the distribution given by Eq. (6), far from γ (n) e,min 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 re-acceleration. In this approximation, the synchrotron and IC peak frequencies increase as ν (n) s ∝ g 2n ν s and roughly ν (n) IC ∝ g n ν IC , respectively. The total emitted synchrotron power per unit volume and frequency at the frequency ν for the power-law electron distribution after n re-accelerations in a turbulent magnetic field is then given following Rybicki & Lightman (2004): with K (n) e the normalisation of the power-law distribution. Keeping the electron number density and magnetic field strength constant, and estimating the peak frequency of the synchrotron component in the δ approximation as ν (n) s [Hz] ≃ 3.7 × 10 6 B [G] γ (n) e,max 2 , the increase in the emitted peak power for the first re-acceleration is roughly 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 re-acceleration. Following Tavecchio, Maraschi, & Ghisellini (1998), it can be written as Article number, page 9 of 18 A&A proofs: manuscript no. zech_lemoine_2021_final The ratio of the Compton dominance between the first reacceleration and the initial acceleration is then given by: (20) The actual spectral distribution of re-accelerated particles deviates from a pure power law, leading to an even more significant effect on the resulting SED. The evolution of a given electron population and its associated broadband emission from an initial shock acceleration to a first and second re-acceleration on similar shocks is shown in Figs. 4 and 5. For this simulation, it was assumed that the number of particles remains constant during the evolution and that radiative cooling can be neglected. The source parameters δ, R, and B were adapted from the solution for 1ES 0229+200 that will be presented below and were left constant for the consecutive accelerations. Under these admittedly simple assumptions, it can be seen that the resulting SED of the electron population is dominated by the last efficient shock acceleration.  Assuming that the emission from 1ES 0229+200 and 1ES 1101-232 is dominated by radiation of a particle population after a first re-acceleration, the SEDs of both sources can be well represented with parameters that are similar to scenario II; in particular a value of γ e,min = 1.8 × 10 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).

Implications for the nature and location of the emission region
The emission region described by the simple one-zone model may represent either a plasma blob moving with relativistic speed through the jet, into which accelerated particles are injected from a leading bow shock (similar to Observed SEDs due to SSC emission from the particle population shown in Fig. 4 after an initial shock acceleration and two re-accelerations. The SEDs correspond to a source at the redshift of 1ES 0229+200 with δ = 50.  Kirk, Rieger, & Mastichiadis 1998), or the region in the relativistic plasma flow downstream from a standing shock.
In both cases, one needs to explain why the emission is restricted mostly to a compact region, resembling the moving or standing VLBI radio knots observed in radio galaxies and blazars. Since in our scenarios synchrotron cooling is very slow, adiabatic cooling or particle escape into regions with even lower magnetic fields may be an explanation. As was discussed above, magnetisation is in fact expected to be amplified in the immediate neighbourhood of the shock. When particles escape from the emission region into the downstream jet, their emission may be shifted to lower fluxes and frequencies. This low-energy emission from the extended jet may be responsible for the radio flux, which we do not include in the SEDs shown here (cf. archival data in Costamante et al. (2018)), and possibly contribute at infrared to ultraviolet frequencies, which may help to better reproduce the SED of 1ES 1218+304 in this range (cf. Fig. 3). Modelling of this extended emission component is beyond the scope of the present study.
The blob slows down in the jet frame once it has swept up a mass M j ∼ E b /(Γ 2 b|j c 2 ) , where E b represents the blob energy (in the jet frame), and Γ b|j the relative blob-jet Lorentz factor. On a timescale t ′ dec (blob frame), the swept-up mass is M j ∼ πR 2 n j m p Γ b|j ct ′ dec , hence the deceleration time can be written t ′ dec ∼ u b / Γ b|j n j m p c 2 R/c, with u b the blob proper energy density. Consequently, t ′ dec is generally larger than t ′ dyn , because the blob ram pressure Γ 2 b|j u b is assumed larger than n j m p c 2 , and Γ b|j is of the order of a few at most. Article number, page 11 of 18 The alternative scenario of acceleration on standing shocks, arising for example from recollimation, may account more naturally for a continuous emission. In the standard picture of an initially Poynting-flux dominated jet, a first recollimation shock is assumed to occur at the transition between the magnetically dominated parabolic base of the jet and a kinetically dominated conical jet at the pc-scale (e.g. Potter 2018). Oblique shocks due to collimation of the relativistic jet by the pressure of the surrounding medium have been considered by, for example, Bromberg & Levinson (2009) to account for the difference between large Doppler factors needed to explain TeV emission and the much smaller values inferred from unification models and super-luminal motion of radio knots. In this framework, acceleration and emission up to VHEs can occur up to several 10 pc from the base of the jet, from very compact emission regions with extensions down to ∼ 10 −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 radio-loud AGN (e.g. Hada et al. 2018, and references therein), and most prominently in the most nearby radio galaxy M 87 at the 'HST-1' radio knot (Asada & Nakamura 2012; Nakamura & Asada 2013), which appears at a de-projected distance of approximately 300 pc from the core, with a radius of the order of 1 pc (Asada & Nakamura 2012).
The typical extensions of the emission regions in our models are of the order of 0.01 pc. To explain the HE emission from extreme-TeV blazars, much larger regions can be ruled out due to the energy requirements, as discussed in the case of 1ES 0229+200 for our scenario II. Besides the above-mentioned predictions from Bromberg & Levinson (2009), an order-ofmagnitude 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.
Re-acceleration on successive standing shocks requires the distance between two shocks to be small compared to the cooling time of particles inside the plasma flow. Following Eq. (2), the most energetic electrons in our models will not significantly cool at distance scales up to about 0.1 Γ j|<s pc (i.e. a few parsecs). This estimate is based on the magnetic field strength in the emission region, which is expected to be amplified compared to regions far from the recollimation shocks; the actual cooling distance may thus be far greater. Formation of internal shocks due to reflections of the converging recollimation shock, proposed by Bromberg & Levinson (2009), may constitute alternative sites of re-acceleration. 2 We note that their definition of σ is slightly different from ours, which should however lead to a result of the same order of magnitude.

Implications for jet power, jet content, and magnetisation
Contrary to other lepto-hadronic models in the literature, which try to explain the observed SEDs as a combination of radiative emission from electrons and protons, emission from protons is negligible in the scenarios proposed here, due to the low magnetic field strengths and small maximum proton Lorentz factors γ p,max . Typical lepto-hadronic models achieve a significant radiative contribution from protons by assuming an acceleration process that is much more efficient for protons than for electrons, leading to u p ≫ u e and γ p,max ≫ γ e,max , combined with high magnetic field strengths, of the order of 1 G to a few 100 G in the case of proton-synchrotron dominated emission. In such models, the jet power is strongly dominated by the contribution from either the magnetic field or the relativistic protons, depending on the dominant hadronic emission mechanism, and is generally close to or even above the Eddington limit (e.g. Mücke & Protheroe 2001;Petropoulou & Mastichiadis 2012;Cerruti et al. 2015, ...).
In our solutions, the magnetic field strength is several orders of magnitude smaller, the number densities of accelerated electrons and protons are equal (n e = n p ), and the jet power is typically a factor 10 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, FSRQ-type blazars show that, at least for those objects, the assumption of a pure electron-proton jet is problematic, since required jet powers are very large compared to the available accretion power and compared to the power that is injected into radio lobes. The addition of an electron-positron pair plasma can lower the required jet power (e.g. Sikora 2016). If this is the case, n e > n p , which would change our constraint on the normalisation of the electron spectrum K e (Eq. 15). The value of γ e,min is expected to decrease by the factor of electron multiplicity. A large multiplicity factor would thus exclude a good representation for the scenarios without substantial reacceleration.
The magnitude of the magnetisation σ of blazar jets is still an open question. When assuming the Blandford-Znajek process for jet launching, this implies an initially strongly Poynting-flux dominated jet, in which magnetic energy is converted to kinetic energy over subparsec to parsec scales. During this conversion, the bulk Lorentz factor increases and the magnetisation drops to unity and below. In the standard picture, this conversion may proceed through differential collimation down to equipartition (i.e. σ ≃ 1), where it becomes inefficient.
Further conversion may be driven by non-standard mechanisms, such as reconnection on magnetohydrodynamic (MHD) instabilities (Komissarov 2011). In the latter case, a magnetisation σ ≪ 1 seems reachable even at small distances from the base of the jet, where most of the radiation is expected to be produced.
In this case, Zdziarski et al. (2015) show that blazar jets can be powered by the Blandford-Znajek mechanism, originate from magnetically arrested disks and can still have magnetisations σ ≪ 1 at the radio core. They also find that in this framework jets are relatively heavy, with their mass dominated by baryons.
Such reconnection processes could lead to particle acceleration. However, if the magnetic field strength in the transition region around σ 1 does not exceed values of the order of 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 multi-zone jet emission model including a realistic description of the evolving reconnection regions, which is well beyond the scope of the current work.
In alternative scenarios where jet launching is described through the Blandford-Payne mechanism, a lower initial magnetisation of the jet may arise more naturally (e.g. Tzeferacos et al. 2009, for simulations in the non-relativistic limit).
Direct measurements of the magnetic field strength B in relativistic jets of blazars and radio galaxies are difficult. They usually rely on measures of the core shift, that is to say, the frequency-dependent position along the jet where synchrotron self-absorption becomes negligible and emission from the VLBI radio core can first be observed. To determine the value of B, usually equipartition between magnetic and particle energy densities is assumed and values of the order of 100 mG are found at a distance of 1 pc from the base of the jet (O'Sullivan & Gabuzda 2010; Kutkin et al. 2014), while Agarwal et al. (2017) find values of the order of 1 mG at the more distant radio cores. Uncertainties and fluctuations in these measurements are however large and deviation from the assumption of equipartition will lead to different values.

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 to 100. These more common objects present at least an order of magnitude lower synchrotron and SSC peak frequencies than the extreme-TeV blazars, which translates into lower values of γ e,max and δ, while the magnetic field strength is usually found to be at least an order of magnitude higher.
In the present framework, the continuous emission from such objects could in principle be accounted for with lower Lorentz factors in the shock normal frame. A higher electron multiplicity would also lead to lower values of γ e,min . Lower values of δ may be due to smaller Lorentz factors in the post-shock region or larger viewing angles of these objects. To reproduce the wider synchrotron and SSC bumps seen in many HBLs, particle spectra require in general a slightly more complex description, usually including a spectral break that cannot be ascribed to radiative cooling of a single electron population alone. These features could possibly be reproduced with a combination of emission from several standing shocks (e.g. a strong first shock) followed by a weaker second shock. In the blob-in-jet scenario, the assumption of a uniform particle density and magnetic field strength in a spherical blob may be too approximate for a realistic description of the particle spectra, hence requiring additional degrees of freedom in the form of an ad hoc spectral break. The acceleration process may also be more complex than what is assumed in our current model. Contrary to extreme TeV blazars, flux and spectral variability are ubiquitous in HBLs, at timescales that can be as short as a few minutes in the observer frame at TeV energies. Additional mechanisms, possibly the interaction of moving shocks with successive standing shocks (e.g. Hervet et al. 2019;Fichet de Clairfontaine et al. 2020, 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.

Conclusions
Ascribing the non-thermal multi-wavelength emission from extreme-TeV blazars to the acceleration of a proton-electron plasma on mildly relativistic shocks introduces additional constraints, arising from the microphysics in the shock region, on the standard one-zone SSC model that is generally applied to describe SEDs of BL Lac objects. A more physically motivated description of the population of radiating electrons can thus be achieved. The co-acceleration of electrons and protons results in a transfer of energy to the electrons, thus shifting their spectral distribution to higher energies. This effect provides a natural explanation for the high values of the minimum electron Lorentz factors required to model the SEDs of extreme-TeV blazars.
A first generic model, allowing for a very high bulk Lorentz factor of the jet plasma upstream of the shock, leads to unique solutions for the set of five extreme-TeV blazar SEDs under study. In a second, more realistic scenario, exploring the two cases of acceleration on recollimation shocks or on the shock wave caused by a plasma blob propagating inside the jet, a good representation is found for most of the SEDs. High minimum Lorentz factors of the electron distribution of γ e,min 10 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 re-acceleration of an initial power-law electron distribution on a second shock can satisfactorily explain the extreme values of the synchrotron and SSC peak frequencies and the hard γ-ray spectrum.
Here we concentrate on the shock crossing conditions and rewrite them in a simple way, in order to derive the post-shock flow velocity and to keep the discussion self-contained. We do so, in particular, by de-boosting to the so-called 'shock normal frame' R n (Begelman & Kirk 1990) in which the flow moves along the shock normal. This means de-boosting by the component of the velocity which is transverse to this shock normal (equivalently, parallel to the shock surface).
We assume that the surface of the recollimation shock makes an angle α with respect to the jet axis. The flow direction, at angle θ j< with respect to this jet axis, thus forms an angle π/2 − θ j< − α with respect to the shock normal. The symbol < is used to distinguish pre-shock from post-shock quantities (indexed with > ). If the shock front is stationary in the source rest frame R s , then the shock normal frame R n moves at velocity β n|s with respect to R s , where β n|s = β j< − β j< · n |s n |s . (B.1) All throughout n |s represents the (unit) direction of the normal to the shock front in the source rest frame. This gives β n|s = β j< cos θ j< − α , hence Γ n|s = 1 − β 2 j< cos 2 θ j< − α −1/2 .
To derive the effective velocity of the jet in this shock normal frame, we note that the component of the jet four-velocity along the shock normal is preserved by the de-boost, hence u j<|n = Γ j< β j< sin θ j< − α . The relative Lorentz factor between the recollimation shock and the incoming plasma in the shock normal frame, which matters for acceleration purposes, can thus be written If Γ j< θ j< − α ≪ 1, the recollimation shock is subrelativistic in the R n frame, otherwise it is relativistic. The postshock downstream velocity can be approximated (in R n ) as β j>|n ≃ β j<n /r in the weakly magnetised limit, with r the compression ratio. Hereafter we assume a mildly relativistic shock with r ≃ 3, implying β j>|n ≃ 0.3, hence Γ j>|n ∼ 1. For reference, the exact shock jump conditions for a strong hydrodynamical shock give β j>|n = 0.29 if Γ j<|n = 3 , β j>|n = 0.26 if Γ j<|n = 2.
Consequently, boosting back to the source frame gives a post-shock plasma four-velocity u j>|s = Γ j>|n β j>|n n |s + Γ j>|n Γ n|s β n|s , (B If Γ j< θ j< − α ≪ 1, the shock is weak, as noted above, and the post-shock plasma flows nearly along the shock surface. In the opposite limit, the recollimation shock is ultra-relativistic in R n (see Eq. (B.2)), the shock is strong and the post-shock flows at sub-or mildly relativistic speeds in R s (see Eq. (B.4)). To describe the intermediate limit, write θ j< − α = κ/Γ j< with κ ∼ O(1). Then Γ j<|n ≃ √ 1 + κ 2 and Γ j>|s ≃ Γ j< / √ 1 + κ 2 . We also note that Eqs. (B.2) and (B.4) give Γ j>|s Γ j<|n = Γ j<|s Γ j>|n ≃ Γ j<|s . From the point of view of modelling, we recall that the post-shock Lorentz factor in the source rest frame (Γ j>|s ) controls the boost factor from the emission region, while the shock Lorentz factor in the shock normal frame (Γ j<|n ) controls the amount of electron heating (i.e. γ e, min ≃ 600Γ j<|n . As the flow passes through the recollimation shock, it is refracted by an angle θ j> − θ j< (in the source rest frame), which can be obtained from the post-shock and pre-shock velocities. In the small-angle approximation, meaning |θ j< |, |θ j> |, |α| ≪ 1, we obtain Therefore, if |θ j< − α| ∼ O 1/Γ j< , the deflection through the recollimation can lead to a substantial change in the Doppler beaming from pre-shock to post-shock values. We ignore this effect in our models, however.

Appendix C: Steady emission at a recollimation shock
As explained in Sikora et al. (1997), the Doppler amplification of the co-moving luminosity takes different forms depending on the geometric situation at hand: a moving blob of matter, or a steady jet. Both are reconciled once the steady jet is interpreted as a train of blobs moving one next to the other. The difference in amplification formula is then related to the actual number of blobs that is seen by an observer at a given time, as a result of relativistic motion. For one moving blob, the general formula is: The subscript |s indicates that the quantities are evaluated in the source rest frame. Henceforth, this subscript is dropped, and quantities written in the co-moving blob or jet frame are primed.
In the case of a steady jet, the standard formula reads 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): where V ′ 1 denotes the co-moving 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 co-moving duration over which a blob is active is δt ′ = d /(Γ b β b ) because δt ′ = δt/Γ b and d = β b δt (δt denoting the emission time in the source rest frame). The apparent width of a blob is δd = r ′ /Γ b , hence (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.