<figure><img src="image_1.png" /></figure>

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

Andreas Zech, Martin Lemoine

► To cite this version:

Andreas Zech, Martin Lemoine. Electron-proton co-acceleration on relativistic shocks in extreme-TeV blazars. Astronomy and Astrophysics - A&A, 2021, 654, pp.A96. 10.1051/0004-6361/202141062. hal-03386490

HAL Id: hal-03386490

https://hal.science/hal-03386490

Submitted on 19 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Electron-proton co-acceleration on relativistic shocks in extreme-TeV blazars

Andreas Zech<fnref n="1" /> and Martin Lemoine<fnref n="2" />

Received 12 April 2021 / Accepted 31 July 2021

## ABSTRACT

**Aims.** 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.

**Methods.** 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.

**Results.** 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.

**Key words.** acceleration of particles – radiation mechanisms: non-thermal – 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

## 1. Introduction

Blazars – active galactic nuclei (AGN) with a jetted relativistic outflow pointing towards us – are commonly classified according to their characteristic double-humped spectral energy distribution (SED; Urry & Padovani 1995; Fossati et al. 1998). At the tip of the so-called blazar sequence (Padovani & Giommi 1995; Ghisellini et al. 2017), extreme blazars form a new population of BL Lacs (Costamante et al. 2001), designated as ‘ultra-high-frequency-peak ed 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 $\gamma$-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 $\gamma$-ray telescope predict about five times more candidate sources that need to be confirmed with future instruments in the TeV range (Costamante 2020). For five of the known extreme-TeV blazars, largely simultaneous SEDs with a good multi-wavelength coverage from the infrared to the TeV band have become available, and they have already provided important clues on the physics of dissipation in these sources (Costamante et al. 2018).

Arriving at a satisfactory description of extreme-TeV sources in the standard one-zone leptonic framework seemingly requires two essential ingredients (e.g., Tavecchio et al. 1998; Ghisellini et al. 2002; Katarzyński et al. 2006; Costamante et al. 2018): (1) lower magnetic fields ($\le 10$ mG) than for more common high-frequency-peak ed 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, $\bar{\gamma}_e \sim 10^3-10^4$. The former requirement avoids an excessive softening of the $\gamma$-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 $\gamma$-ray spectrum and high peak frequency. Alternative models have been devised (e.g., Aharonian 2000; Böttcher et al. 2008; Lefa et al. 2011; Tavecchio 2014; Cerruti et al. 2015; Chhotray et al. 2017; Tavecchio & Sobacchi 2020), but they all come at the price of greater complexity and most often invoke other extreme parameter values, which often lack experimental or physical justification. For instance, lepto-hadronic models can provide reasonable SEDs with less extreme values of $\bar{\gamma}_e$, but at the expense of a high jet power, mostly carried by the accompanying population of ultra-relativistic protons or the magnetic field (e.g., Cerruti et al. 2015).

In the present study we take at face value the low magnetisation inferred from the SEDs of extreme blazars and

<footnote n="1">¹ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France
e-mail: andreas.zech@obspm.fr</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, CNRS – Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France
e-mail: lemoine@iap.fr</footnote>
propose revisiting the one-zone model, adopting well-motivated microphysical prescriptions for particle acceleration, under the assumption that electrons are co-accelerated with protons in relativistic internal or recollimation shocks. Two key motivations of our study are the following: (1) at low magnetisation, relativistic shock acceleration is expected to dominate the physics of dissipation and of the generation of non-thermal power laws, and (2) electrons can be efficiently preheated in the transition layer of electron-ion relativistic shocks, up to a fraction of the thermal energy of shocked ions, implying large minimum Lorentz factors for the electron population, as seemingly required by phenomenology (e.g., Vanthieghem et al. 2020, and references therein). We further aim at making our approach as self-contained as possible, by fixing the microphysical parameters describing the acceleration process from first principle studies of relativistic shock models. Our main objective is to provide a critical discussion of such models to reproduce the SEDs of extreme blazars, in a modern context. This contrasts with generic one-zone models, which instead parameterise the accelerated particle population and adjust the parameters to the observed SEDs.

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

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

### 2.1. A large $\bar{\gamma}_e$

We define the Lorentz factor $\bar{\gamma}_e$ as the average over the population,

$$ \bar{\gamma}_e = \frac{1}{n_e} \int d\gamma_e \gamma_e \frac{dn_e}{d\gamma_e}, \qquad (1) $$

with $n_e = \int d\gamma_e dn_e/d\gamma_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 $\gamma_{e,b}$. For such a distribution of Lorentz factors, $\bar{\gamma}_e \sim \gamma_{e,min} (\gamma_{e,b}/\gamma_{e,min})^{2-s_1}$ if $s_1 < 2 < s_2$, and $\bar{\gamma}_e \sim \gamma_{e,min}$ if $2 < s_1 < s_2$. The energy density $u_e$ of the population can then be written $u_e \sim \bar{\gamma}_e n_e m_e c^2$.

Therefore, to convert a population of electrons (or pairs) of initial Lorentz factor $\gamma_0$ to the above power-law spectrum necessitates the dissipation of an energy reservoir that is $\bar{\gamma}_e/\gamma_0$ times larger than the initial particle energy. In this sense, the large inferred value of $\bar{\gamma}_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), $\bar{\gamma}_e \sim 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 $\gamma_e = 10^4 \gamma_{e4}$ cool through synchrotron radiation on a length scale

$$ l_c \approx 8\Gamma_j B_{-2}^{-2} \gamma_{e4}^{-1} \text{pc}, \qquad (2) $$

as measured in the source rest frame, and expressed in terms of $\Gamma_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 $\Gamma_j^4 n_e R^3 u_B \gamma_{e,min}^{\gamma_{e,min}} \gamma_{e,max}^{\gamma_{e,max}}$ for a single power law of index $s < 3$ extending from $\gamma_{e,min}$ to $\gamma_{e,max}$, where $u_B$ denotes the magnetic energy density, $n_e$ the electron density, $R$ the jet radius, and $\Gamma_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 $\gamma_{sh}$ (and velocity $\beta_{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} \equiv \gamma_{sh}\beta_{sh} \ge 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 \approx 0.19 u_{sh}^2$, while for a strong, relativistic pair shock, $k_B T_e/m_e c^2 \approx 0.24 u_{sh}$. For reference, $k_B T_e/m_e c^2 \approx 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 = \gamma_{sh} (\gamma_{sh} - 1) n_p m_p c^2$ ($n_p$ is understood as a proper density here), corresponding to an energy-per-particle $\langle E_p \rangle = (\gamma_{sh} - 1) m_p c^2$, and similarly for the electrons with $p \leftrightarrow e$. Efficient preheating of the electrons means that $\langle E_e \rangle$ becomes a substantial fraction of $\langle E_p \rangle$, instead of $\langle E_e \rangle / \langle E_p \rangle = m_e/m_p \ll 1$, which would be expected if both species were to satisfy the fluid shock crossing conditions independently from each other.
The mechanism(s) through which electrons gain energy at the expense of ions in shock fronts remain somewhat elusive. In the relativistic, weakly magnetised limit, one promising scenario is that in which electrons undergo collisionless Joule heating, driven by a longitudinal electric field and by the effective gravity that results from the slowdown of the microturbulence that is self-generated in the shock precursor (see e.g., Lemoine et al. 2019c; Vanthieghem et al., in prep.). Phenomenological models of gamma-ray burst afterglows provide a nice empirical confirmation of this energy transfer between ions and electrons, as the electron energy fraction is almost always found to be within a factor of a few from the ion energy fraction (e.g., Kumar & Zhang 2015, and references therein).

In the absence of a definitive model, we can rely on particle-in-cell (PIC) simulations to provide estimates for the average electron Lorentz factor. Generally speaking, the preheating appears more efficient at lower magnetisations $\sigma$ (Sironi et al. 2013), likely because of a larger precursor and because of the presence of intense Weibel-like microturbulence. For reference, we define $\sigma$ 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 $\gamma_{sh} \ge 3-5$ – than in the mildly relativistic regime $\gamma_{sh} \sim 1-3$. For example, Crumley et al. (2019) measure $T_e \approx 0.3T_p$ and $k_B T_p \approx 0.2 m_p c^2$ for a mildly relativistic sub-luminal shock of velocity $\beta_{sh} = 0.8$, corresponding to $\gamma_{sh} = 1.7$, at magnetisation $\sigma \sim 0.007$ (here defined downstream). This electronic temperature corresponds to a mean Lorentz factor $\bar{\gamma}_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 \sim 0.1 \gamma_{sh} m_p c^2$, or a mean Lorentz factor $\bar{\gamma}_e \sim 600 \gamma_{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 $\epsilon_e$ of the energy that is incoming into the shock, for example $\epsilon_e \approx 0.001$ for the case with $\beta_{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 $\bar{\gamma}_e$. In the relativistic limit, $\epsilon_e \ge 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 $\bar{\gamma}_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 $\bar{\gamma}_e$, but also a hard spectrum with index $s \sim 2$. In the following, we thus adopt $\epsilon_e = 0.1$ and $\bar{\gamma}_e \sim 600 \gamma_{sh}$ as fiducial values, corresponding to relativistic, weakly magnetised electron-ion shocks. We also assume one electron per ion, as larger electron multiplicities are expected to lead to lower mean energy per particle, although this case has not been explicitly probed by numerical experiments.

## 2.3. Fixing parameters from the microphysics of shock acceleration

In the fully relativistic regime, meaning a shock Lorentz factor $\gamma_{sh} \ge 3-5$, particle acceleration appears restricted to the weakly magnetised limit, $\sigma \le 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 $\gamma_{sh}$) of the perpendicular magnetic field components in the shock rest frame. In the mildly relativistic limit, however, such effects disappear and sub-luminal configurations become about as likely. Consequently, particle acceleration can become efficient, even at relatively high magnetisations (see e.g., Crumley et al. 2019).

In the regime of low magnetisation that we are interested in, the size of the precursor is governed by the scattering of particles in the electromagnetic micro-turbulence that they themselves excite through beam-plasma instabilities (e.g., Plotnikov et al. 2018), and which is responsible for the shock transition itself. This general problem has been studied in some detail in the limit of negligible initial magnetisation (e.g., Lemoine et al. 2019a,b,c; Pelletier et al. 2019). The microturbulence is mostly excited by the current filamentation instability (often termed Weibel), which generates an intense magnetic field on skin depth scales $\lambda_p$, with $\lambda_p \equiv c (4\pi n_p e^2/m_p)^{-1/2} \sim 10^7 n_p^{1/2} \text{ cm} (n_{p0} = n_p/1 \text{ cm}^{-3})$. The effective magnetisation carried by this turbulence, written $\epsilon_B$, generally saturates at values $\epsilon_B \sim 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 $\sigma_{\text{rad}}$, is expected to be smaller than $10^{-2}$, but at least as large as the initial $\sigma$ value.

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

At moderate magnetisation, $10^{-5} \lesssim \sigma \lesssim 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 $\epsilon_B \sim 0.01$ and a length scale $\sim O(\lambda_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 \approx 2-2.3$. For definiteness, we assume $s \approx 2.2$ in our application to extreme-TeV blazars, corresponding to mildly relativistic shocks. The minimum Lorentz factor is $\gamma_{e,\min} \approx \bar{\gamma}_e \approx 600 \gamma_{sh}$ as discussed before and the maximal Lorentz factor $\gamma_{e,\max}$ is determined by the conjunction of a number of constraints: (1) the age limit, $t_{\text{acc}} \lesssim d/(\Gamma_j \beta_j c)$, where $t_{\text{acc}}$ represents the acceleration timescale to $\gamma_{e,\max}$ and $d$ the distance from the jet base, $\Gamma_j$ (resp. $\beta_j$) the
jet Lorentz factor (resp. velocity); (2) the radiative loss limit, $t_{acc} \lesssim t_{syn}$, written here as a competition with the synchrotron loss timescale; (3) the scattering limit (for super-luminal shocks) $t_{scatt} \lesssim 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} \lesssim r$, is comparable to the age limit, if $r \sim \Theta_j d$ and the jet opening angle $\Theta_j \sim 1/\Gamma_j$.

The acceleration timescale is given by $t_{acc} \approx t_{scatt}$ in relativistic shocks, where $t_{scatt}$ denotes the scattering timescale, which departs from the Bohm regime ($t_{scatt} \propto r_L$) at weak magnetisation (Vanthieghem et al. 2020, and references therein). Consequently, the radiative limit can be written $\gamma_{e,max} \sim 10^7 (n_e/1 \text{ cm}^{-3})^{-1/6}$, which does not depend on the strength of the magnetic field in the acceleration region, contrary to the usual relation $\gamma_{e,max} \propto B^{-1/2}$ obtained for Bohm scaling.

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

$$ \gamma_{e,max} \lesssim \frac{\gamma_{e,min}}{\sqrt{\sigma}}. \qquad (3) $$

Its main effect is to reduce the dynamic range over which particles can be accelerated. Its dependence on $\sigma$ derives from the scalings $t_{scatt} \propto (\gamma_e/\gamma_{sh})^2$ and $r_{L,0} \propto (\gamma_e/\gamma_{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 $\epsilon_B = 10^{-2}$, $\sigma = 10^{-5}$, $\gamma_{sh} = 3$, $n_e = 1 \text{ cm}^{-3}$. Then, the scattering constraint gives $\gamma_{e,max} \sim 2 \times 10^5$, while the radiative and age constraints both lead to $\gamma_{e,max} \sim 10^7$. The severity of the scattering constraint results from the relatively “low” value of $\gamma_{e,min}$, here 1800; below, we discuss models in which $\gamma_{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, $\gamma_{p,min} \sim \gamma_{sh}$, a maximal Lorentz factor determined by the conjunction of the age and scattering constraints above and the same spectral index as for the electrons.

## 2.4. Shock origin and geometry

### 2.4.1. Shocks of normal incidence

The interaction of a jet with an obstacle generally triggers a double shock configuration, including a forward and a reverse shock.

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

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

Alternatively, if an ‘obstacle’ overtakes the jet, as in the ‘blob-in-jet’ scenario, the Doppler factor will be at least as large as that of the jet itself. If the blob ram pressure exceeds that of the jet, the Lorentz factor of the emission region is given by that of the blob, $\Gamma_b$, while the shock Lorentz factor $\gamma_{sh}$ is set as before by the relative Lorentz factor between the blob and the jet. From $\gamma_{sh} = \Gamma_j\Gamma_b(1 - \beta_j\beta_b)$, $\Gamma_j \gg 1$, $\Gamma_b \gg 1$ and $\Gamma_b \gg \Gamma_j$, it follows

$$ \gamma_{sh} \sim \Gamma_b/(2\Gamma_j). \qquad (4) $$

Mildly relativistic shocks can be attained for sufficiently large ratios of $\Gamma_b$ to $\Gamma_j$. We consider such a configuration in Sect. 3.2.

### 2.4.2. Oblique shocks

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

In brief, these constraints remain applicable to oblique shocks once they are expressed in the so-called shock normal rest frame, in which the (un-shocked) flow impinges on the shock at normal incidence (Begelman & Kirk 1990). The Lorentz factor $\gamma_{sh}$, which corresponds to the Lorentz factor of the shock with respect to the upstream plasma in this shock normal frame, is written $\Gamma_{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 $\Gamma_{j<}$, through Eq. (B.2): $\gamma_{sh} = \Gamma_{j<}/\Gamma_{n|s}$, where $\Gamma_{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 $|\theta_{j<} - \alpha|$ between the flow incidence and the shock surface ($\theta_{j<}$ and $\alpha$, 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/\Gamma_j$, one expects $|\theta_{j<} - \alpha| \sim 1/\Gamma_j$ in a simple re-confinement scenario (Nalewajko & Sikora 2009), and we therefore write $|\theta_{j<} - \alpha| \sim \kappa/\Gamma_{j<}$. Then, we obtain $\gamma_{sh} \approx \sqrt{1 + \kappa^2}$. The shock is thus mildly relativistic for $\kappa$ of the order of unity. The angle $|\theta_{j<} - \alpha|$ 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 $\Gamma_{j>} \approx \Gamma_{n|s} \approx \Gamma_{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 $|\theta_{j<} - \alpha|$ 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 $\delta^4$ (cf. Sikora et al. 1997), with $\delta$ the bulk Doppler factor. For the sake of completeness, we provide more details on this issue in Appendix C.

### 2.4.3. Interaction with multiple shocks

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

As a particle crosses a shock front, its energy in the local frame of rest of the plasma is increased on average by the relative Lorentz factor between pre-shock and post-shock flows. Average is understood here over the orientation of the particle momentum in the pre-shock flow. For an oblique shock, the relative Lorentz factor $\Gamma_{rel}$ can be written in terms of normal shock frame quantities $\Gamma_{rel} \approx 0.7\Gamma_{j<\ln}$, where the prefactor holds for mildly relativistic shocks with $\Gamma_{j<\ln} \ge 2 - 3$. In terms of source frame quantities, $\Gamma_{rel} \approx 0.5\Gamma_{j<}/\Gamma_{j>}$ (in analogy with Eq. (4)).

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

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

$$G(\gamma_{>}, \gamma_{<}) = \frac{s-1}{g\gamma_{<}} \left(\frac{\gamma_{>}}{g\gamma_{<}}\right)^{-s} \Theta(\gamma_{>} - g\gamma_{<}), \quad (5)$$

which provides the probability distribution of the outgoing particle Lorentz factor $\gamma_{>} \approx \Gamma_{rel}^{2/3}$ as discussed above. The post-shock particle distribution can be obtained through the convolution of the pre-shock distribution with this propagator, which gives, after n shock crossings (White 1985; Achterberg 1990; Schneider 1993)

$$\frac{dN_{>}^{(n)}}{d\gamma_{>}} = \frac{(s-1)^{n+1}}{n! g^n \gamma_{\min}} \left(\frac{\gamma_{>}}{g^n \gamma_{\min}}\right)^{-s} \ln \left(\frac{\gamma_{>}}{g^n \gamma_{\min}}\right)^n. \quad (6)$$

This has two important consequences: (1) the spectrum hardens because of reacceleration; (2) the effective injection Lorentz factor, which was $\gamma_{\min}$ at the first shock, has become $g^n\gamma_{\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\gamma$, which provides a reasonable estimate for $\bar{\gamma}_e$ is found at $\gamma = g^n\gamma_{\min} \exp[n/(s-1)]$, which indicates a rapid rise of the effective $\bar{\gamma}_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, $\gamma_{e,\min}$ is now given by $g^n \gamma_{e,\min}$, with $\gamma_{e,\min} \approx 600\gamma_{sh}$ and the spectrum is given by Eq. (6) above.

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

The standard one-zone SSC model has usually nine free parameters: the strength of the uniform magnetic field B, the size of the spherical emission region R, its Doppler factor $\delta$, the minimum, maximum and break Lorentz factors of the electron distribution $\gamma_{e,\min}$, $\gamma_{e,\max}$, $\gamma_{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 $\gamma = 1^1$. Even with a well-sampled multi-wavelength dataset, this model remains highly degenerate.

Tavecchio et al. (1998) have studied the constraints on these parameters for a set of six observables extracted from the SED: the peak frequencies and luminosities of the synchrotron and inverse Compton (IC) components, $\nu_{syn}$, $\nu_{IC}$, $\nu_{syn}L_{syn}(\nu_{syn})$, $\nu_{IC}L_{IC}(\nu_{IC})$, and the indices in the differential photon flux $F_\gamma$ before and after the synchrotron peak, $\alpha_1$ and $\alpha_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 $\gamma_{e,\min}$, $\gamma_{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

<footnote n="1">¹ All model parameters refer to a frame co-moving with the emission region, i.e. the plasma ‘blob’.</footnote>
by considering a power-law electron distribution with an exponential cutoff of the form

$$ \frac{dN_e}{d\gamma} = K_e \gamma^{-s} e^{-\gamma/\gamma_{e,max}}. \qquad (7) $$

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

For an initially chosen value of $\delta$, first the maximum electron Lorentz factor $\gamma_{e,max}$ can be determined from the IC peak frequency $\nu_{IC}$. In the Klein-Nishina regime, which is relevant for the datasets we want to investigate in the following section, this relation is given following Tavecchio et al. (1998):

$$ \gamma_{e,max} \approx \frac{1}{g(\alpha_1, \alpha_2)} \frac{1+z}{\delta} \frac{hv_{IC}}{m_e c^2}, \qquad (8) $$

with

$$ g(\alpha_1, \alpha_2) = \exp\left(\frac{1}{\alpha_1 - 1} + \frac{1}{2(\alpha_2 - \alpha_1)}\right). \qquad (9) $$

Setting $s \approx 2.2$ then implies $\alpha_1 \approx 0.6$ as $\alpha = (s - 1)/2$ because the particles do not effectively cool through synchrotron radiation. The value of $\alpha_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 $\gamma_{e,max}$ known, the magnetic field strength $B$ in the radiation region can be estimated from the value of $\nu_{syn}$, giving

$$ B = 27 \, \text{mG} \frac{1+z}{\delta} \frac{\nu_{syn}}{10^{17} \, \text{Hz}} \left(\frac{\gamma_{e,max}}{10^6}\right)^{-2}. \qquad (10) $$

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

Provided the maximum electron energy is set by the scattering constraint, following Eq. (3), and the minimum electron energy can be estimated, at least broadly, from the low-frequency part of the SED and the spectral slope in the GeV range, the magnetisation of the plasma flow can be estimated as $\sigma \approx (\gamma_{e,min}/\gamma_{e,max})^2$.

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

$$ R = \left\{ 2 \frac{f(\alpha_1, \alpha_2)}{\delta^4} \frac{[\nu_{syn} L_{v_{syn}}]^2}{\nu_{IC} L_{v_{IC}} B^2 c} \left( \frac{3\delta m_e c^2}{4\gamma_{e,max} h\nu_{syn}} \right)^{(1-\alpha_1)} \right\}^{1/2}, \qquad (11) $$

with

$$ f(\alpha_1, \alpha_2) = \frac{1}{1 - \alpha_1} + \frac{1}{\alpha_2 - 1}. \qquad (12) $$

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

$$ \gamma_c \approx 7.7 \times 10^8 \left( \frac{B}{1 \, \text{mG}} \right)^{-2} \left( \frac{R}{3 \times 10^{16} \, \text{cm}} \right)^{-1} \qquad (13) $$

versus $\gamma_{e,max}$. Given the low magnetisation required for efficient shock acceleration, the assumption holds in our application. IC

cooling is not significant either in the absence of strong Compton dominance.

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

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

$$ u_p = \int_{\gamma_{p,min}}^{\gamma_{p,max}} d\gamma K_p \gamma^{-s+1} e^{-\gamma/\gamma_{p,max}} m_p c^2 = \frac{u_B}{\sigma_{rad}}, \qquad (14) $$

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

$$ K_e = K_p \frac{\int_{\gamma_{p,min}}^{\gamma_{p,max}} d\gamma \gamma^{-s} e^{-\gamma/\gamma_{p,max}}}{\int_{\gamma_{e,min}}^{\gamma_{e,max}} d\gamma \gamma^{-s} e^{-\gamma/\gamma_{e,max}}}. \qquad (15) $$

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

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

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

$$ L_j \approx 2\pi R^2 \beta c \Gamma^2 (u_B + u_e + u_p + u_\gamma), \qquad (16) $$

where $u_\gamma$ represents the energy density contained in the SSC radiation field.

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

### 3. Application to extreme-TeV blazars

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

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

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

The co-acceleration model was declined into three scenarios to find the physically most meaningful representations of the selected SEDs. In a first scenario, we study a generic version with a large value of $\gamma_{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 $\gamma_{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 1: Generic shock model with large $\gamma_{e,min}$

We thus consider here a model with a large value of $\gamma_{e,min}$. In light of our discussion in Sect. 2, such a large value can

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

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., ~$\delta$), means that the jet encounters an obstacle moving at small velocity in the source rest frame (see Sect. 2.4.1). However, this requires the jet ram pressure to exceed that of the obstacle in the source frame, which is precisely opposite to the condition that allows the obstacle to penetrate the jet (e.g., Araudo et al. 2010; Barkov et al. 2010; Bosch-Ramon et al. 2012). Consequently, none of the above possibilities currently receives clear justification. For this reason, we regard this first scenario as a case study. To fix the parameters, we assume that $\gamma_{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 $\delta$, 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 $\gamma_{e,min}$, of a few $10^4$, in agreement with previous modelling attempts that can be found in the literature, where $\gamma_{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. $\sigma_{rad} = \sigma$). 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 $\delta$ 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 $\delta = 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 $\delta$ 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 $\gamma_{e,min}$ by an order of magnitude, in addition to adjustments of $\delta$ and $R$ (cf. Fig. 2).

Allowing for these adjustments, satisfactory solutions can be found for all five SEDs. The full application of the constraints
<figure><img src="image_3.png" /><figcaption>Fig. 2. SSC model for the SED of 1ES 1218+304 in scenario I.</figcaption></figure>

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 $\delta = 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 $\sigma \ll 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.

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

In this second scenario, we consider the more physically grounded case of electron preheating and acceleration at internal (Sect. 2.4.1) or recollimation shocks (Sect. 2.4.2). In both cases, far smaller values of the shock Lorentz factor in the shock normal frame $\Gamma_{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 $\Gamma_b$ traversing the jet, the resulting shock Lorentz factor $\gamma_{sh} \equiv \Gamma_{j<|n}$ is given by Eq. (4). Assuming a Doppler factor of $\delta$ for the observed emission from the blob, a mildly relativistic shock with $\Gamma_{j<|n} \approx 3$ would result for a jet Lorentz factor $\Gamma_j \approx \delta/12$, when assuming a closely aligned jet, corresponding to $\Gamma_b \approx \delta/2$, such that $\Gamma_j$ would be of the order of a few. Such a combination of a highly relativistic emission component responsible for HE emission and a mildly relativistic jet is frequently proposed to explain the discrepancy in observed bulk Lorentz factors in the radio and HE bands (Ghisellini et al. 2005; Henri & Pelletier 1991; Sol et al. 1989). In this case, shock Lorentz factors much higher than a few can only be achieved for unrealistically high blob Lorentz factors or for non-relativistic jets.

Considering the second case of acceleration on recollimation shocks, since the bulk Lorentz factor of the downstream, radiating plasma is $\Gamma_{j>|s} \ge \delta/2$, the shock Lorentz factor must verify $\Gamma_{j<|n} \approx \Gamma_{j<|s}/\Gamma_{j>|s} \le 2\Gamma_{j<|s}/\delta$ (cf. Appendix B). If the bulk Lorentz factor of the plasma flow upstream of the recollimation shock $\Gamma_{j<|s}$ is of a similar strength as the typical bulk Doppler factors derived from SED modelling, the shock Lorentz factor

$\Gamma_{j<|n}$ should be close to unity. The default value of $\Gamma_{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 $\gamma_{e,min} \approx 600 \Gamma_{j<|n} \approx 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 $\delta \approx 50$, as for the solutions in scenario I, the decrease in $\gamma_{e,min}$ for a fixed $\gamma_{e,max}$ leads to a quadratic decrease in the inferred upstream magnetisation $\sigma$ (cf. Eq. (3)). Since the particle density is inversely proportional to the magnetisation downstream of the shock (cf. Eq. (14)), the condition $\sigma_{rad} = \sigma$ 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 $\sigma \le \sigma_{rad} \le 0.01$, with 0.01 representing the typical effective magnetisation at the shock, $\sigma_{rad}$ the magnetisation in the radiation region and $\sigma$ the pre-shock magnetisation. We also stress that the values of $\sigma$ that are provided in Table 1 are only indicative, in the sense that they are derived from the respective values of $\gamma_{e,min}$ and $\gamma_{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 $\sigma_{rad}$. Consequently, the true physical magnetisation to be regarded is $\sigma_{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 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 $\sigma_{rad}$ introduces a degeneracy in the sense that similar solutions can be found by increasing $\delta$ and decreasing $\sigma_{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 $\gamma$-ray data points, as a result of the smaller value of $\gamma_{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 $\gamma_{e,min}$ remains fixed at its default value, the only way around this problem is to increase $\delta$ to a value well above 100, in order to shift the low-energy break to sufficiently high frequencies.

Otherwise, satisfactory solutions for 1ES 0229+200 and 1ES 1101–232 can be obtained only if the shock Lorentz factor $\Gamma_{j<|n}$ is increased to large values, thus recovering a scenario
<table><caption>**Table 1.** Parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.</caption><thead><tr><th></th><th>1ES 0229+200</th><th>1ES 0347-121</th><th>RGB J0710+591</th><th>1ES 1101-232</th><th>1ES 1218+304</th></tr></thead><tbody><tr><td colspan="6"><strong>Scenario I</strong></td></tr><tr><td>δ</td><td>50</td><td>30</td><td>40</td><td>50</td><td>60</td></tr><tr><td>R<sub>src</sub> [10<sup>16</sup>cm]</td><td>1.8</td><td>4.7</td><td>2.0</td><td>3.0</td><td>1.1</td></tr><tr><td>B [mG]</td><td>4.4</td><td>3.7</td><td>7.9</td><td>8.9</td><td>8.4</td></tr><tr><td>σ<sub>rad</sub> (≡ σ)</td><td>1.0 × 10<sup>-4</sup></td><td>1.4 × 10<sup>-4</sup></td><td>5.0 × 10<sup>-4</sup></td><td>1.6 × 10<sup>-3</sup></td><td>9.8 × 10<sup>-5</sup></td></tr><tr><td>γ<sub>e,min</sub></td><td>3.0 × 10<sup>4</sup></td><td>1.8 × 10<sup>4</sup></td><td>2.4 × 10<sup>4</sup></td><td>3.0 × 10<sup>4</sup></td><td>5.0 × 10<sup>3</sup></td></tr><tr><td>γ<sub>e,max</sub></td><td>2.9 × 10<sup>6</sup></td><td>1.5 × 10<sup>6</sup></td><td>1.1 × 10<sup>6</sup></td><td>7.6 × 10<sup>5</sup></td><td>5.0 × 10<sup>5</sup></td></tr><tr><td>n<sub>e</sub> [cm<sup>-3</sup>]</td><td>2.9 × 10<sup>-2</sup></td><td>2.7 × 10<sup>-2</sup></td><td>2.8 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-2</sup></td><td>0.69</td></tr><tr><td colspan="6"><strong>Scenario II</strong></td></tr><tr><td>δ</td><td>50</td><td>50</td><td>50</td><td>50</td><td>60</td></tr><tr><td>R<sub>src</sub> [10<sup>16</sup>cm]</td><td>2.3</td><td>1.6</td><td>1.9</td><td>4.5</td><td>1.4</td></tr><tr><td>B [mG]</td><td>4.4</td><td>6.2</td><td>9.8</td><td>8.9</td><td>8.4</td></tr><tr><td>σ<sub>rad</sub></td><td>1.9 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-4</sup></td><td>1.2 × 10<sup>-3</sup></td><td>4.0 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>σ<sub>rad</sub>/σ</td><td>11</td><td>29</td><td>270</td><td>64</td><td>11</td></tr><tr><td>γ<sub>e,min</sub></td><td>1.2 × 10<sup>4</sup></td><td>1.8 × 10<sup>3</sup></td><td>1.8 × 10<sup>3</sup></td><td>6.0 × 10<sup>3</sup></td><td>1.8 × 10<sup>3</sup></td></tr><tr><td>γ<sub>e,max</sub></td><td>2.9 × 10<sup>6</sup></td><td>9.1 × 10<sup>5</sup></td><td>8.6 × 10<sup>5</sup></td><td>7.6 × 10<sup>5</sup></td><td>5.0 × 10<sup>5</sup></td></tr><tr><td>n<sub>e</sub> [cm<sup>-3</sup>]</td><td>4.2 × 10<sup>-2</sup></td><td>0.83</td><td>0.24</td><td>2.1 × 10<sup>-2</sup></td><td>1.3</td></tr><tr><td colspan="6"><strong>Scenario III</strong></td></tr><tr><td>δ</td><td>50</td><td>50</td><td>—</td><td>50</td><td>—</td></tr><tr><td>R<sub>src</sub> [10<sup>16</sup>cm]</td><td>1.3</td><td>0.9</td><td>—</td><td>7.0</td><td>—</td></tr><tr><td>B [mG]</td><td>4.4</td><td>6.2</td><td>—</td><td>4.5</td><td>—</td></tr><tr><td>σ<sub>rad</sub></td><td>1.0 × 10<sup>-4</sup></td><td>9.0 × 10<sup>-5</sup></td><td>—</td><td>4.0 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>σ<sub>rad</sub>/σ</td><td>270</td><td>23</td><td>—</td><td>720</td><td>—</td></tr><tr><td>γ<sub>e,min</sub></td><td>1.8 × 10<sup>3</sup></td><td>1.8 × 10<sup>3</sup></td><td>—</td><td>1.8 × 10<sup>3</sup></td><td>—</td></tr><tr><td>γ<sub>e,max</sub></td><td>2.9 × 10<sup>6</sup></td><td>9.1 × 10<sup>5</sup></td><td>—</td><td>7.6 × 10<sup>5</sup></td><td>—</td></tr><tr><td>n<sub>e</sub> [cm<sup>-3</sup>]</td><td>0.13</td><td>0.29</td><td>—</td><td>8.2 × 10<sup>-3</sup></td><td>—</td></tr></tbody></table>

Notes. The observables used as input for the modelling and parameters of the proton population can be found in Table D.1.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** SSC model for the SED of 1ES 1218+304 in scenario II.</figcaption></figure>

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 $\Gamma_{j<|s|} \approx 20$ is necessary, while in the case of 1ES 1101-232, $\Gamma_{j<|s|} \approx 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 $\Gamma_{j<|s|} \approx 10^3$ for 1ES 0229+200, while for 1ES 1101-232, it would result in

$\Gamma_{j<|s|} \approx 500$. Such large jet Lorentz factors are problematic given that the usually assumed values for less extreme blazars, based on observations with very-long-baseline interferometry (VLBI) and on models of relativistic jets, are of the order of a few to a few tens (Pushkarev et al. 2009; Saikia et al. 2016).

For acceleration on internal shocks caused by blob-jet interaction, the required values of $\Gamma_{j<|s|}$ 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 $\delta \lesssim \Gamma_b$).

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

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

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

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

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

$$P_{\text{syn}}^{(n)}(v) \propto K_{e}^{(n)} \frac{1}{s^{(n)} + 1} \Gamma\left(\frac{s^{(n)}}{4} + \frac{19}{12}\right) \Gamma\left(\frac{s^{(n)}}{4} - \frac{1}{12}\right) \times \left(\frac{2\pi v m_e c}{3eB}\right)^{-\frac{s^{(n)}-1}{2}}, \quad (17)$$

with $K_e^{(n)}$ the normalisation of the power-law distribution.

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

$$\frac{v_s^{(1)} P_{\text{syn}}^{(1)}(v_s^{(1)})}{v_s^{(0)} P_{\text{syn}}^{(0)}(v_s^{(0)})} \approx 7.1 \left(\frac{\gamma_{e,max}}{10^6}\right)^{0.1}. \quad (18)$$

The numerical factor in the above equation depends weakly on the value of $\gamma_{e,min}$, which we set again to $1.8 \times 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 et al. (1998), it can be written as

$$\kappa^{(n)} \equiv \left(\frac{L_{\text{IC}}}{L_{\text{syn}}}\right)^{(n)} \propto f(\alpha_1^{(n)}, \alpha_2) v_s^{(n)} P_{\text{syn}}^{(n)}(v_s^{(n)}) \times \left(\frac{3 m_e c^2}{4 \times 3.7 \times 10^6 h}\right)^{(3-s^{(n)})/2} \left(\frac{\delta}{B \gamma_{e,\text{max}}^{(n)}\delta}\right)^{(3-s^{(n)})/2} . \quad (19)$$

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

$$\frac{\kappa^{(1)}}{\kappa^{(0)}} \approx 2.6 \left( \frac{\delta}{50} \frac{1 \text{ mG}}{B} \frac{10^6}{\gamma_{e,\text{max}}} \right)^{0.05}. \quad (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

<figure><img src="image_5.png" /><figcaption>Fig. 4. Particle spectra for the same electron and proton population after the initial shock acceleration and after a first and second re-acceleration on successive shocks.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Observed SEDs due to SSC emission from the particle population shown in Fig. 4 after an initial shock acceleration and two re-accelerations. The SEDs correspond to a source at the redshift of 1ES 0229+200 with $\delta = 50$.</figcaption></figure>

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 $\delta$, $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 $\gamma_{e,min} = 1.8 \times 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 $v_{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).
<figure><img src="image_7.png" /><figcaption>Fig. 6. SSC model for the SED of 1ES 0229+200 in scenario III.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. SSC model for the SED of 1ES 1101–232 in scenario III.</figcaption></figure>

# 4. Discussion

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

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

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

The blob-in-jet scenario can accommodate reasonably well a baryon loaded medium with low magnetisation if the blob

interacts with the jet at large distances from the jet base. However, the blob emission timescale is constrained both by the expansion of the blob and its slowdown, due to its interaction with the jet. We assume that the expansion of the blob takes place at a lateral velocity $v_⊥$ in the source frame, with $v_⊥/c \sim O(\Theta_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'_{dyn} \sim R/(v_⊥\Gamma_b)$, corresponding to an observer timescale of $\Delta t_{obs} = t'_{dyn}/\delta \sim R/(v_⊥\Gamma_b\delta)$. If $v_⊥/c \sim \Theta_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 $\Theta_j \ll 1/\Gamma$ 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_j \sim E_b/(\Gamma_{blj}^2 c^2)$, where $E_b$ represents the blob energy (in the jet frame), and $\Gamma_{blj}$ the relative blob-jet Lorentz factor. On a timescale $t'_{dec}$ (blob frame), the swept-up mass is $M_j \sim \pi R^2 n_j m_p \Gamma_{blj} c t'_{dec}$, hence the deceleration time can be written $t'_{dec} \sim [u_b / (\Gamma_{blj} 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 $\Gamma_{blj}^2 u_b$ is assumed larger than $n_j m_p c^2$, and $\Gamma_{blj}$ is of the order of a few at most.

The alternative scenario of acceleration on standing shocks, arising for example from recollimation, may account more naturally for a continuous emission. In the standard picture of an initially Poynting-flux dominated jet, a first recollimation shock is assumed to occur at the transition between the magnetically dominated parabolic base of the jet and a kinetically dominated conical jet at the pc-scale (e.g., Potter 2018). Oblique shocks due to collimation of the relativistic jet by the pressure of the surrounding medium have been considered by, for example, Bromberg & Levinson (2009) to account for the difference between large Doppler factors needed to explain TeV emission and the much smaller values inferred from unification models and super-luminal motion of radio knots. In this framework, acceleration and emission up to VHEs can occur up to several 10 pc from the base of the jet, from very compact emission regions with extensions down to $\sim 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 M87 at the ‘HST-1’ radio knot (Asada & Nakamura 2012; Nakamura & Asada 2013), which appears at a de-projected distance of approximately 300 pc from the core, with a radius of the order of 1 pc (Asada & Nakamura 2012).

The typical extensions of the emission regions in our models are of the order of 0.01 pc. To explain the HE emission from extreme-TeV blazars, much larger regions can be ruled out due to the energy requirements, as discussed in the case of 1ES 0229+200 for our scenario II. Besides the above-mentioned predictions from Bromberg & Levinson (2009), an order-of-magnitude estimate of the location of the emission regions in our models can be made by comparing their radii with the assumed jet opening angle $\Theta_j$, since $\Theta_j \sim R/d$, with $d$ the distance of the emission region from the base of the jet.
Typical values for $\Theta_j$, based on VLBI data from the MOJAVE programme for γ-loud blazars are found to follow roughly $\Theta_j \simeq 0.3/\Gamma_{jl<s}$ (Pushkarev et al. 2009), while Zdziarski et al. (2015) find $\Theta_j \simeq 0.1/\Gamma_j$ in a different sample. The latter note however that the opening angle is related to the magnetisation as $\Theta_j \simeq s\sigma^{1/2}/\Gamma_{jl<s}$ with $s \lesssim 1$, which would suggest a relation closer to $\Theta_j \simeq 0.01/\Gamma_{jl<s}$ for the typical values of $\sigma$ found in our scenarios<fnref n="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 \Gamma_{jl<s}$ pc (i.e. a few parsecs). This estimate is based on the magnetic field strength in the emission region, which is expected to be amplified compared to regions far from the recollimation shocks; the actual cooling distance may thus be far greater. Formation of internal shocks due to reflections of the converging recollimation shock, proposed by Bromberg & Levinson (2009), may constitute alternative sites of re-acceleration.

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

Contrary to other lepto-hadronic models in the literature, which try to explain the observed SEDs as a combination of radiative emission from electrons and protons, emission from protons is negligible in the scenarios proposed here, due to the low magnetic field strengths and small maximum proton Lorentz factors $\gamma_{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 \gg u_e$ and $\gamma_{p,max} \gg \gamma_{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 $\sigma \ll 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 $\gamma_{e,min}$ is expected to decrease by the factor of electron multiplicity. A large multiplicity factor would thus exclude a good representation for the scenarios without substantial re-acceleration.

The magnitude of the magnetisation $\sigma$ 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. $\sigma \simeq 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 $\sigma \ll 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 $\sigma \ll 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 $\sigma \lesssim 1$ does not exceed values of the order of 1 G, the emission from those particles would contribute mainly to the extended radio and optical emission of the jet, not to the high-frequency part, because the average electron Lorentz factor decreases with decreasing $\sigma$ and the produced spectra become quite soft for $\sigma \lesssim 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} \lesssim 1$ G, $\sigma_{trans} \lesssim 1$ and a steep particle spectrum with $\gamma_{e,min} \sim 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} \sim 50 (R_{trans}/10^{16} \text{ cm})^3 (B_{trans}/1 \text{ G})^4 (1/\sigma_{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.

<footnote n="2">² We note that their definition of $\sigma$ is slightly different from ours, which should however lead to a result of the same order of magnitude.</footnote>
Uncertainties and fluctuations in these measurements are how-
ever 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 γ<sub>e,min</sub>, typically of the order of 1–100. These more common objects present at least an order of magnitude lower synchrotron and SSC peak frequencies than the extreme-TeV blazars, which translates into lower values of γ<sub>e,max</sub> 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 multi-
plicity would also lead to lower values of γ<sub>e,min</sub>. 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, par-
ticle spectra require in general a slightly more complex descrip-
tion, 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 realis-
tic 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 vari-
ability are ubiquitous in HBLs, at timescales that can be as
short as a few minutes in the observer frame at TeV ener-
gies. Additional mechanisms, possibly the interaction of mov-
ing 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. Alterna-
tively, depending on the magnetisation in the emission region,
scenarios based on particle acceleration through reconnection or
turbulence may provide a suitable description for the variable
emission from such objects.

5. Conclusions

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

A first generic model, allowing for a very high bulk Lorentz factor of the jet plasma upstream of the shock, leads to unique solutions for the set of five extreme-TeV blazar SEDs under study. In a second, more realistic scenario, exploring the two cases of acceleration on recollimation shocks or on the shock wave caused by a plasma blob propagating inside the jet, a good

representation is found for most of the SEDs. High minimum
Lorentz factors of the electron distribution of γ<sub>e,min</sub> ≥ 10³ arise
as a natural consequence of the assumed microphysics.

For the SEDs of the two objects with the hardest γ-ray
spectra, 1ES 0229+200 and 1ES 1101-232, good representa-
tions are found for acceleration on moving shocks, but also for
the acceleration on recollimation shocks when assuming that the
particle distribution that dominates the emission results from
re-acceleration on a second shock. In the latter case, the spec-
tral 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 satisfac-
torily explain the extreme values of the synchrotron and SSC
peak frequencies and the hard γ-ray spectrum.

***Acknowledgements.*** 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

Agarwal, A., Mohan, P., Gupta, A. C., et al. 2017, MNRAS, 469, 813

Aharonian, F. A. 2000, New Astron., 5, 377

Araudo, A. T., Bosch-Ramon, V., & Romero, G. E. 2010, A&A, 522, A97

Asada, K., & Nakamura, M. 2012, ApJ, 745, L28

Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80

Barkov, M. V., Aharonian, F. A., & Bosch-Ramon, V. 2010, ApJ, 724, 1517

Begelman, M. C., & Kirk, J. G. 1990, ApJ, 353, 66

Biteau, J., Prandini, E., Costamante, L., et al. 2020, Nat. Astron., 4, 124

Bosch-Ramon, V., Perucho, M., & Barkov, M. V. 2012, A&A, 539, A69

Böttcher, M., Dermer, C. D., & Finke, J. D. 2008, ApJ, 679, L9

Bromberg, O., & Levinson, A. 2009, ApJ, 699, 1274

Cerruti, M., Boisson, C., & Zech, A. 2013, A&A, 558, A47

Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910

Chhotray, A., Nappo, F., Ghisellini, G., et al. 2017, MNRAS, 466, 3544

Costamante, L. 2020, MNRAS, 491, 2771

Costamante, L., Ghisellini, G., Giommi, P., et al. 2001, A&A, 371, 512

Costamante, L., Bonnoli, G., Tavecchio, F., et al. 2018, MNRAS, 477, 4257

Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, MNRAS, 485, 5105

Dmytriiev, A., Sol, H., & Zech, A. 2021, MNRAS, 505, 2712

Dubal, M. R., & Pantano, O. 1993, MNRAS, 261, 203

Duffell, P. C., & MacFadyen, A. I. 2014, ApJ, 791, L1

Fichet de Clairfontaine, G., Meliani, Z., Zech, A., & Hervet, O. 2021, A&A, 647, A77

Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433

Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837

Ghisellini, G., Celotti, A., & Costamante, L. 2002, A&A, 386, 833

Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401

Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011, MNRAS, 414, 2674

Ghisellini, G., Righi, C., Costamante, L., & Tavecchio, F. 2017, MNRAS, 469, 255

Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ, 860, 141

Henri, G., & Pelletier, G. 1991, ApJ, 383, L7

Hervet, O., Williams, D. A., Falcone, A. D., & Kaur, A. 2019, ApJ, 877, 26

H.E.S.S. Collaboration (Abdalla, H., et al.) 2020, Nature, 582, 356

Inoue, T., Asano, K., & Ioka, K. 2011, ApJ, 734, 77

Katarzyński, K., Ghisellini, G., Tavecchio, F., Gracia, J., & Maraschi, L. 2006, MNRAS, 368, L52

Kawazura, Y., Schekochihin, A. A., Barnes, M., et al. 2020, PhRvX, 10, 041050

Kirk, J. G., & Duffy, P. 1999, JPhG, 25, R163

Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452

Komissarov, S. S. 2011, MmSAI, 82, 95

Komissarov, S. S., & Falle, S. A. E. G. 1997, MNRAS, 288, 833

Kumar, P., & Zhang, B. 2015, PhR, 561, 1

Kutkin, A. M., Sokolovsky, K. V., Lisakov, M. M., et al. 2014, MNRAS, 437, 3396

Lefa, E., Rieger, F. M., & Aharonian, F. 2011, ApJ, 740, 64

Lemoine, M. 2013, MNRAS, 428, 845

Lemoine, M., & Pelletier, G. 2010, MNRAS, 402, 321

Lemoine, M., Li, Z., & Wang, X.-Y. 2013, MNRAS, 435, 3009

Lemoine, M., Pelletier, G., Gremillet, L., & Plotnikov, I. 2014, EL, 106, 55001
Lemoine, M., Vanthieghem, A., Pelletier, G., & Gremillet, L. 2019a, PhRvE, 100, 033209

Lemoine, M., Pelletier, G., Vanthieghem, A., & Gremillet, L. 2019b, PhRvE, 100, 033210

Lemoine, M., Gremillet, L., Pelletier, G., & Vanthieghem, A. 2019c, PhRvL, 123, 035101

Martí, J. M., Perucho, M., & Gómez, J. L. 2016, ApJ, 831, 163

Meli, A., & Biermann, P. L. 2013, A&A, 556, A88

Mücke, A., & Protheroe, R. J. 2001, APh, 15, 121

Nakamura, M., & Asada, K. 2013, ApJ, 775, 118

Nalewajko, K., & Sikora, M. 2009, MNRAS, 392, 1205

O'Sullivan, S. P., & Gabuzda, D. C. 2010, ASPC, 427, 207

Padovani, P., & Giommi, P. 1995, ApJ, 444, 567

Pelletier, G., Lemoine, M., & Marcowith, A. 2009, MNRAS, 393, 587

Pelletier, G., Gremillet, L., Vanthieghem, A., & Lemoine, M. 2019, PhRvE, 100, 013205

Perucho, M. 2013, EPJWC, 61, 02002

Petropoulou, M., & Mastichiadis, A. 2012, MNRAS, 426, 462

Petropoulou, M., Sironi, L., Spitkovsky, A., & Giannios, D. 2019, ApJ, 880, 37

Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238

Pope, M. H., & Melrose, D. B. 1994, PASAu, 11, 175

Potter, W. J. 2018, MNRAS, 473, 4107

Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2009, A&A, 507, L33

Rybicki, G. B., & Lightman, A. P. 2004, Radiative Process in Astrophysics (Hoboken: Cfa Harvard, Wiley-VCH)

Saikia, P., Körding, E., & Falcke, H. 2016, MNRAS, 461, 297

Santana, R., Barniol, Duran R., & Kumar, P. 2014, ApJ, 785, 29

Sari, R., & Piran, T. 1995, ApJ, 455, L143

Sbarrato, T., Foschini, L., Ghisellini, G., & Tavecchio, F. 2011, AdSpR, 48, 998

Schneider, P. 1993, A&A, 278, 315

Sikora, M. 2016, Galax, 4, 12

Sikora, M., Madejski, G., Moderski, R., & Poutanen, J. 1997, ApJ, 484, 108

Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103

Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75

Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54

Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411

Tavecchio, F. 2014, MNRAS, 438, 3255

Tavecchio, F., & Sobacchi, E. 2020, MNRAS, 491, 2198

Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608

Tzeferacos, P., Ferrari, A., Mignone, A., et al. 2009, MNRAS, 400, 820

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Vanthieghem, A., Lemoine, M., Plotnikov, I., et al. 2020, Galax, 8, 33

Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840

White, R. L. 1985, ApJ, 289, 698

Zdziarski, A. A., Sikora, M., Pjanka, P., & Tchekhovskoy, A. 2015, MNRAS, 451, 927

Zech, A., Cerruti, M., & Mazin, D. 2017, A & A, 602, A25

Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, PhRvL, 122
# Appendix A: Electron heating in reconnection layers or magnetised turbulence

To produce electron spectra with large $\bar{\gamma}_e$, reconnection or turbulent acceleration must proceed in the relativistic regime, because they need to dissipate an energy reservoir ~ $\bar{\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 \bar{\gamma}_e/\gamma_0 \gg 1$, clearly at odds with the apparent requirement of a sub-equipartition magnetic field in the radiation region. As mentioned earlier, one may contemplate the possibility that such values are reached at the core of the jet in a region of high magnetisation, provided the radiation produced in this region does not dominate. Such dissipation generally leads, however, to an approximate equipartition between particles and magnetic fields, leading to a magnetisation $\sigma \sim O(1)$, which remains too large to account for the spectral properties of extreme blazars, whence the need for some additional mechanism to bring down this magnetisation to the desired sub-equipartition level by the time the electron power law is shaped.

In an electron-proton plasma, the initial energy density is carried by the protons, hence reconnection can proceed in the relativistic regime from the point of view of electrons ($\sigma_e > 1$, but in the sub-relativistic regime from the point of view of ions ($\sigma_p \le 1$). Here, $\sigma_e$ (respectively, $\sigma_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 $\sigma$ 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 $\sigma \ge 1$, suggest $\bar{\gamma}_e \sim 1 + 0.03\sigma^{3/2}m_p/m_e$ for an initially cold electron population, hence values $\sigma \ge 10$ appear needed to reach $\bar{\gamma}_e \sim 10^3$. This result remains essentially unchanged if the electrons are initially hot ($T_e \gg m_e$). At lower magnetisations, the electron spectra take a more complex form, composed by part of their initial Maxwellian distribution at low Lorentz factors and a non-thermal population emerging out of a secondary Maxwellian at large Lorentz factor $\gamma_e^{\text{hi}}$. From the results of Ball et al. (2018) and Werner et al. (2018), we infer $\gamma_e^{\text{hi}} \sim (0.1-0.3)\sigma m_p/m_e$. Since $\bar{\gamma}_e \lesssim \gamma_e^{\text{hi}}$ for such electron spectra, we recover the above result: in other words, $\sigma \gg 1$ is required to obtain $\bar{\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 \sim 0.25+0.25\sigma^{1/2}/(10+\sigma)^{1/2}$ in the regime $\sigma > 0.01$, which, when combined with $\bar{\gamma}_p \approx 1+0.5\sigma$ at $\sigma \le 1$ (Ball et al. 2018), gives $\bar{\gamma}_e \approx 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 \approx 4$ for $\sigma \le 1$ and becoming as hard as $s \sim 2$ for $\sigma \gg 1$. All this clearly indicates that reconnection does not fulfil the requirements for a successful model for extreme blazars under our present assumptions.

The partition of energy between electrons and ions in MHD turbulence is a rather convoluted problem, which depends on a number of parameters, not only the initial magnetisation and species temperatures, but also the nature (compressive or Alfvénic) of the turbulence itself (e.g. Kawazura et al. 2020). As of today, there exist only a few kinetic simulations of ion-electron turbulence in the regime of interest and we rely on the results of Zhdankin et al. (2019). These authors have performed a scan of simulations which span a broad range of magnetisations, from $\sigma_p \sim 10^{-3}$ (hence, $\sigma_e \sim 1$) up to the fully relativistic

regime $\sigma_p \ge 1$ ($\sigma_e \sim 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 \approx m_e$, the electrons acquire about 10% of the ion energy gain. If $\gamma_p - 1 \sim \sigma$, as expected for near-complete transfer of the magnetic energy into particles, this means $\bar{\gamma}_e \approx 0.1\sigma m_p/m_e$, a result comparable to what is observed in reconnection layers. Hence, $\bar{\gamma}_e \ge 10^3$ can be achieved for large magnetisations only.

# Appendix B: Jump conditions at a recollimation shock

Here we concentrate on the shock crossing conditions and rewrite them in a simple way, in order to derive the post-shock flow velocity and to keep the discussion self-contained. We do so, in particular, by de-boosting to the so-called ‘shock normal frame’ $\mathcal{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 $\alpha$ with respect to the jet axis. The flow direction, at angle $\theta_{j<}$ with respect to this jet axis, thus forms an angle $\pi/2 - (\theta_{j<} - \alpha)$ 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 $\mathcal{R}_s$, then the shock normal frame $\mathcal{R}_n$ moves at velocity $\beta_{n|s}$ with respect to $\mathcal{R}_s$, where

$$ \beta_{n|s} = \beta_{j<} - (\beta_{j<} \cdot n_{|s}) n_{|s}. \quad (B.1) $$

All throughout $n_{|s}$ represents the (unit) direction of the normal to the shock front in the source rest frame. This gives $\beta_{n|s} = \beta_{j<} \cos(\theta_{j<} - \alpha)$, hence $\Gamma_{n|s} = [1 - \beta_{j<}^2 \cos^2(\theta_{j<} - \alpha)]^{-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} = \Gamma_{j<}\beta_{j<} \sin(\theta_{j<} - \alpha)$. 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

$$ \Gamma_{j<|n} = \sqrt{1 + \Gamma_{j<}^2 \beta_{j<}^2 \sin^2(\theta_{j<} - \alpha)} = \frac{\Gamma_{j<}}{\Gamma_{n|s}}. \quad (B.2) $$

If $|\theta_{j<} - \alpha| \ll 1$, the recollimation shock is sub-relativistic in the $\mathcal{R}_n$ frame, otherwise it is relativistic. The post-shock downstream velocity can be approximated (in $\mathcal{R}_n$) as $\beta_{j>|n} \approx \beta_{j<|n}/r$ in the weakly magnetised limit, with $r$ the compression ratio. Hereafter we assume a mildly relativistic shock with $r \approx 3$, implying $\beta_{j>|n} \approx 0.3$, hence $\Gamma_{j>|n} \sim 1$. For reference, the exact shock jump conditions for a strong hydrodynamical shock give $\beta_{j>|n} = 0.29$ if $\Gamma_{j<|n} = 3$, $\beta_{j>|n} = 0.26$ if $\Gamma_{j<|n} = 2$.

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

$$ u_{j>|n} = \Gamma_{j>|n} \beta_{j>|n} n_{|s} + \Gamma_{j>|n} \Gamma_{n|s} \beta_{n|s}, \quad (B.3) $$

with corresponding Lorentz factor

$$ \Gamma_{j>|n} = \Gamma_{j>|n} \Gamma_{n|s} \approx \frac{\Gamma_{j<}}{\sqrt{1 + \beta_{j<}^2 \Gamma_{j<} \sin^2(\theta_{j<} - \alpha)}}. \quad (B.4) $$
If $\Gamma_{j<} |\theta_{j<} - \alpha| \ll 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 $\mathcal{R}_n$ (see Eq. (B.2)), the shock is strong and the post-shock flows at sub- or mildly relativistic speeds in $\mathcal{R}_s$ (see Eq. (B.4)). To describe the intermediate limit, write $|\theta_{j<} - \alpha| = \kappa/\Gamma_{j<}$ with $\kappa \sim \mathcal{O}(1)$. Then $\Gamma_{j<|n} \approx \sqrt{1 + \kappa^2}$ and $\Gamma_{j>|s} \approx \Gamma_{j<}/\sqrt{1 + \kappa^2}$.

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

As the flow passes through the recollimation shock, it is refracted by an angle $\theta_{j>} - \theta_{j<}$ (in the source rest frame), which can be obtained from the post-shock and pre-shock velocities. In the small-angle approximation, meaning $|\theta_{j<}|, |\theta_{j>}|, |\alpha| \ll 1$, we obtain

$$|\theta_{j>} - \theta_{j<}| \approx \frac{|\theta_{j<} - \alpha|}{\sqrt{3}}. \qquad (B.5)$$

Therefore, if $|\theta_{j<} - \alpha| \sim \mathcal{O}(1/\Gamma_{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:

$$F_v = \frac{1+z}{d_L^2} \left. \frac{dE}{dt_{\text{obs}}d\nu d\Omega} \right|_s = \frac{1+z}{d_L^2} \delta^3 \int dV' j_{v'}' \qquad (C.1)$$

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

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

$$F_v = \frac{1+z}{d_L^2} \frac{\delta^2}{\Gamma_b} \int dV' j_{v'}' \qquad (C.2)$$

because $dt_{\text{obs}} = dt$, versus $dt_{\text{obs}} = dt(1 - \beta_b \mu)$ for the moving blob. However, describing the steady jet as a train of blobs leads to Sikora et al. (1997):

$$F_v = \frac{1+z}{d_L^2} \delta^3 N_{\text{b,eff}} \int_{V'_1} dV' j_{v'}', \qquad (C.3)$$

where $V'_1$ denotes the co-moving volume of one blob, of co-moving size $r'$. The total volume of the emission region (in the

source frame) is assumed to be composed of $N_b$ blobs at a given time. Because the emission time of a given blob is compressed by $(1 - \beta_b \mu) = 1/(\Gamma_b \delta)$ in the observer direction, the effective number of blobs that the observer sees at any given time is not $N_b$, but only $N_{\text{b,eff}} = N_b/ (\Gamma_b \delta)$. Then both descriptions indeed match one another, given that $V = V'/\Gamma_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_\parallel$ (the length of the train of blobs in the source rest frame). The co-moving duration over which a blob is active is $\delta t' = d_\parallel/(\Gamma_b \beta_b)$ because $\delta t' = \delta t/\Gamma_b$ and $d_\parallel = \beta_b \delta t$ ($\delta t$ denoting the emission time in the source rest frame). The apparent width of a blob is $\delta d_\parallel = r'/\Gamma_b$, hence

$$N_b = \frac{d_\parallel}{\delta d_\parallel} \approx \Gamma_b^2 \beta_b \frac{\delta t'}{r'}, \qquad (C.4)$$

and consequently $N_{\text{b,eff}} \sim \delta t'/r'$. If $\delta t' \sim r'$, as in the usual blob model, then $N_{\text{b,eff}} \sim 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_\parallel$ can be fixed by other considerations than the assumption $\delta t' \sim 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_\parallel$, 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_\parallel \sim r'/\Theta_j$ and $\Theta_j \sim 1/\Gamma_b$, one recovers $N_{\text{b,eff}} \sim 1$.

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

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

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

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.
<table><caption>**Table D.1.** Observables and proton parameters used for the modelling of the extreme blazar datasets in scenarios I, II, and III.</caption><thead><tr><th></th><th>1ES 0229+200</th><th>1ES 0347-121</th><th>RGB J0710+591</th><th>1ES 1101-232</th><th>1ES 1218+304</th></tr></thead><tbody><tr><td>z</td><td>0.140</td><td>0.188</td><td>0.125</td><td>0.186</td><td>0.184</td></tr><tr><td>L<sub>Edd</sub> [10<sup>46</sup> erg s<sup>-1</sup>]</td><td>18</td><td>1.3</td><td>2.2</td><td>13</td><td>1.4</td></tr><tr><td>t<sub>var,obs</sub> [s]</td><td>~ 10<sup>7</sup></td><td>-</td><td>-</td><td>-</td><td>~ 10<sup>5</sup></td></tr><tr><td>ν<sub>syn</sub> [Hz]</td><td>3 × 10<sup>18</sup></td><td>4 × 10<sup>17</sup></td><td>6 × 10<sup>17</sup></td><td>2 – 4 × 10<sup>17</sup></td><td>2 × 10<sup>17</sup></td></tr><tr><td>ν<sub>IC</sub> [Hz]</td><td>1 × 10<sup>27</sup></td><td>3 × 10<sup>26</sup></td><td>3 × 10<sup>26</sup></td><td>2.5 × 10<sup>26</sup></td><td>2 × 10<sup>26</sup></td></tr><tr><td>νF<sub>v</sub>(ν<sub>syn</sub>) [erg cm<sup>-2</sup> s<sup>-1</sup>]</td><td>2.0 × 10<sup>-11</sup></td><td>5.0 × 10<sup>-12</sup></td><td>1.3 × 10<sup>-11</sup></td><td>2.0 × 10<sup>-11</sup></td><td>1.0 × 10<sup>-11</sup></td></tr><tr><td>νF<sub>v</sub>(ν<sub>IC</sub>) [erg cm<sup>-2</sup> s<sup>-1</sup>]</td><td>5.0 × 10<sup>-12</sup></td><td>5.0 × 10<sup>-12</sup></td><td>6.3 × 10<sup>-12</sup></td><td>6.3 × 10<sup>-12</sup></td><td>1.5 × 10<sup>-11</sup></td></tr><tr><td colspan="6">scenario I</td></tr><tr><td>γ<sub>p,min</sub></td><td>50</td><td>30</td><td>40</td><td>50</td><td>8</td></tr><tr><td>γ<sub>p,max</sub></td><td>4.9 × 10<sup>3</sup></td><td>2.5 × 10<sup>3</sup></td><td>1.8 × 10<sup>3</sup></td><td>1.3 × 10<sup>3</sup></td><td>840</td></tr><tr><td>L<sub>jet</sub> [10<sup>44</sup> erg s<sup>-1</sup>]</td><td>3.6</td><td>5.1</td><td>2.1</td><td>2.9</td><td>8.3</td></tr><tr><td colspan="6">scenario II</td></tr><tr><td>γ<sub>p,min</sub></td><td>20</td><td>3</td><td>3</td><td>10</td><td>3</td></tr><tr><td>γ<sub>p,max</sub></td><td>1.5 × 10<sup>3</sup></td><td>280</td><td>87</td><td>160</td><td>250</td></tr><tr><td>L<sub>jet</sub> [10<sup>44</sup> erg s<sup>-1</sup>]</td><td>3.6</td><td>5.7</td><td>2.1</td><td>2.9</td><td>9.2</td></tr><tr><td colspan="6">scenario III</td></tr><tr><td>γ<sub>p,min</sub></td><td>3.0</td><td>3.0</td><td>-</td><td>3.0</td><td>-</td></tr><tr><td>γ<sub>p,max</sub></td><td>300.</td><td>320.</td><td>-</td><td>47.</td><td>-</td></tr><tr><td>L<sub>jet</sub> [10<sup>44</sup> erg s<sup>-1</sup>]</td><td>3.7</td><td>2.6</td><td>-</td><td>2.5</td><td>-</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. D.1.** SSC model for the SED of RGB J0710+591 in scenario I.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. D.3.** SSC model for the SED of 1ES 1101-232 in scenario I.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. D.2.** SSC model for the SED of 1ES 0347-121 in scenario I.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. D.4.** SSC model for the SED of RGB J0710+591 in scenario II.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. D.5. SSC model for the SED of 1ES 0347-121 in scenario II.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. D.7. SSC model for the SED of 1ES 0229+200 in scenario II.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. D.6. SSC model for the SED of 1ES 1101-232 in scenario II.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. D.8. SSC model for the SED of 1ES 0347-121 in scenario III.</figcaption></figure>