# Extreme TeV BL Lacs: a self-consistent stochastic acceleration model

Alberto Sciaccaluga1,2, Fabrizio Tavecchio1
1INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46, I–23807 Merate,
Italy
2Dipartimento di Fisica, Università degli studi di Genova, Via Dodecaneso 33,
16146, Italy
E–mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Lately a specific kind of blazars drew the attention of the gamma-ray
astronomy community: the extreme TeV BL Lacs, blazars that present an
extremely energetic and hard emission at very high-energy. Explaining their
features is still an open challenge, in fact the most used phenomenological
models have difficulties to satisfactorily reproduce their SED. Based on a
scenario we have recently proposed, we suppose that the non-thermal particles
are firstly accelerated by a jet recollimation shock, which induces turbulence
in the rest of the jet. Non-thermal particles are further accelerated by the
turbulence, which hardens the particle spectra and accordingly the radiative
emission. Given the physical properties of the plasma, as inferred by emission
models, we expect a strong impact of the accelerating particles on the
turbulence. Assuming isotropy and homogeneity, the interaction between non-
thermal particles and turbulence and their spectra is modeled solving a system
of two non-linear, coupled Fokker-Planck equations, while the radiative
emission is calculated through the Synchrotron Self Compton model. The
emission predicted by our model is then compared with the prototype extreme
TeV BL Lac object 1ES 0229+200 and the parameters obtained to reproduce its
SED are in line with the expectations.

###### keywords:

galaxies: jets — radiation mechanisms: non-thermal — shock waves — gamma–rays:
galaxies

##  1 Introduction

Blazars are active galactic nuclei with a relativistic jet (produced in the
supermassive black hole vicinity and propagating far outside the host galaxy)
pointing toward the Earth (Romero et al., 2017; Blandford et al., 2019).
Thanks to the relativistic beaming their emission is dominated by the jet non-
thermal component, characterized by an emission ranging from radio to gamma-
ray. The spectral energy distribution displays two broad peaks (Ghisellini et
al., 2017), attributed to synchrotron and inverse Compton mechanisms (but the
high-energy peak is also interpreted in terms of hadronic processes, see e.g.
Böttcher et al. 2013). The maximum of the second component shifts from the MeV
band, for the most powerful flat spectrum radio quasars, to the TeV band for
BL Lacs. Generally, the variability of the blazar emission is extreme, both in
amplitude and timescale (Aharonian et al. 2006, Bonnoli et al. 2011).

Among blazars, extreme highly peaked BL Lac objects (EHBLs; Costamante et al.
2001; Biteau et al. 2020) are a peculiar class which recently emerged because
of their unique features. In particular we focus on the extreme TeV BL Lacs, a
subclass of the EHBLs which displays extreme properties: i) the second SED
peak beyond $1\ $\mathrm{TeV}$$; ii) a hard sub-TeV intrinsic spectrum; iii)
in most cases, the TeV emission is stable over years. These characteristics
are hardly reproducible by a phenomenological model based only on diffusive
shock acceleration by a single shock. In fact the sub-TeV spectrum hardness
requires the particle spectrum to have spectral index $p<2$ (supposing
$dN/d\gamma\propto\gamma^{-p}$), but with a single shock acceleration $p\geq
2$. In order to obtain harder spectra, other acceleration scenarios have been
proposed, based on a hard electron distribution with a large minimum energy
(Tavecchio et al., 2009), a maxwellian-like electron distribution (Lefa et
al., 2011), a beam of high-energy hadrons (Essey & Kusenko, 2010), internal
absorption (Aharonian et al., 2008) or emission from a large-scale jet
(Böttcher et al., 2008).

Leptonic models which leave the acceleration process unspecified show that the
mean magnetic field and, consequently, the magnetization of extreme TeV Bl
Lacs are low (Biteau et al., 2020), therefore the diffusive shock acceleration
should work efficiently (Sironi et al., 2015). Recently Zech & Lemoine (2021)
proposed a model based on diffusive shock acceleration by multiple shocks.
After an initial shock due to the recollimation of the jet by the external
gas, several shocks could form, as confirmed by the 2D numerical simulations
on jet recollimation (Fichet de Clairfontaine et al., 2021). Non-thermal
particles are then accelerated multiple times and they can reproduce the hard
spectra observed for the extreme TeV Bl Lacs. However, recent 3D simulations
of the jet recollimation (Gourgouliatos & Komissarov, 2018; Matsumoto et al.,
2021) show that if the magnetization is low, after the first shock the jet
becomes unstable and turbulent, preventing the formation of additional shocks.
Based on this new evidence we have proposed a different scenario: the non-
thermal particles are firstly accelerated by the recollimation shock, then
they are further accelerated by the turbulence formed in the downstream
region. This hybrid scenario (see also Kundu et al. 2021) based on diffusive
shock acceleration and stochastic acceleration for extreme TeV BL Lacs has
been invoked in Tavecchio et al. 2022, but there are some important caveats,
in particular regarding the damping of the turbulence, that can modify the
expected emission. Indeed, since the energy density of the non-thermal
particles is much more larger than the magnetic energy density, a strong
damping of the turbulence is expected. In these conditions the standard
Kolmogorov spectrum adopted in Tavecchio et al. 2022 is not self-consistent.
We note that this is feature common to models for blazar emission based on
stochastic acceleration, for which strong damping is generally expected but
not self-consistently determined (see for example Kakuwa et al. 2015).

In this paper we therefore extend the simple treatment of Tavecchio et al.
2022, introducing a self-consistent model capable of describing the interplay
between the turbulence damping and the non-thermal particles acceleration.

The paper is structured as follows: the phenomenological model based on
turbulence is introduced in Sect. 2, then the details of the numerical method
are presented (Sect. 3) and compared with the existing data (Sect. 4). In
Sect. 5 we report our conclusions.

Throughout the paper, the following cosmological parameters are assumed:
$H_{0}=70{\rm\;km\;s}^{-1}{\rm\;Mpc}^{-1}$, $\Omega_{\rm M}=0.3$,
$\Omega_{\Lambda}=0.7$.

##  2 The model

Along the lines of the scenario proposed in Tavecchio et al. 2022, the model
presented in this paper is based on a double acceleration mechanism, an hybrid
of diffusive shock acceleration and stochastic acceleration due to turbulence
(see also Kundu et al. 2020).

We suppose that the non-thermal particles are first accelerated by a
recollimation shock, formed in the jet when the external pressure is greater
than the jet pressure. As already explained, the shock acceleration process is
a deterministic process and the output spectra of the non-thermal particles is
a power law, with spectral index $p\geq 2$ (Blandford & Eichler, 1987; Sironi
& Spitkovsky, 2011).

The non-thermal particles injected by the recollimation shock can be further
energized interacting with the turbulence, hardening the distribution and
therefore the emitted spectrum. Non-thermal particles are supposed to be
accelerated mostly through gyroresonant interaction with the weakly turbulent
Alfvén waves propagating parallel or antiparallel to the mean magnetic field,
therefore the phase speed is equal to the Alfvén speed (for non-resonant
acceleration see Bresci et al. 2022). This acceleration process is stochastic
and in order to study the equilibrium spectra of both the non-thermal particle
spectra and the turbulence it is necessary to consider the time evolution of
their phase-space densities through two momentum diffusion equations.

Assuming isotropy and homogeneity, the momentum diffusion equation can be
highly simplified, obtaining the classical Fokker-Planck equation. Therefore
our goal is to solve a system of two coupled non-linear Fokker-Planck
equations, which describe both the interaction and the spectra of the non-
thermal particles and of turbulence. The emission recorded by the observer is
produced in the shock downstream region and comprises non-thermal particles at
different stages of acceleration, from the shock to the top of the emitting
region, which is modelled as a cylinder of radius $R$ and length $10\,R$. All
quantities appearing in equations below are expressed in the downstream jet
frame, we will move to the observer frame only for the quantities linked to
the radiative emission.

The first equation describes the time evolution of $f(p)$, the momentum
distribution of the non-thermal particles, from here on supposed to be
relativistic electrons, since we are considering a leptonic model:

$\frac{\partial f}{\partial t}=\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}D_{p}\frac{\partial f}{\partial p}+p^{2}\left(\frac{\partial
p}{\partial t}\right)_{\text{rad}}f\right]+\frac{f}{t_{\text{esc}}}+I_{f}$ (1)

where $p$ is the particle momentum, $D_{p}$ the momentum diffusion
coefficient, $\left(\partial p/\partial t\right)_{\text{rad}}$ the cooling
coefficient due to radiative emission and $I_{f}$ is the particle injection
rate. Starting from momentum distribution, it is possible to calculate the
electron energy density, in fact $n(\gamma)=4\pi p^{2}f(p)m_{e}c$, where
$m_{e}$ is the electron mass. Supposing only parallel and antiparallel
propagating MHD waves interacting with electrons, the acceleration time can be
estimated as (Miller & Roberts, 1995; Kakuwa, 2016):

$t_{\text{acc}}=\left[\frac{2\beta_{w}^{2}c}{U_{B}r_{g}}\int_{k_{\text{res}}}\frac{W_{B}(k)}{k}dk\right]^{-1}$
(2)

where $\beta_{w}=v_{a}/c$ is the wave velocity in unit of c,
$U_{B}=B^{2}/8\pi$ is the magnetic energy density, $r_{g}=\gamma
m_{e}c^{2}/eB$ is the Larmor radius, $W_{B}(k)\approx W(k)/2$ is the magnetic
component of the energy density of turbulence fields per unit of $k$
(indicated with $W(k)$) and $k_{\text{res}}=2\pi/r_{g}$ is the resonant
wavenumber. Note that an electron with a fixed $\gamma$ can interact only with
Alfvén waves whose wavelength is smaller than its Larmor radius. From the time
acceleration the momentum diffusion coefficient can be obtained using
$D_{p}=p^{2}/2t_{\text{acc}}$. Neglecting inverse Compton and adiabatic
cooling for simplicity, the cooling time is equal to:

$t_{\text{cool}}=\frac{6\pi m_{e}c}{\sigma_{T}B^{2}\gamma}$ (3)

where $B$ is the mean magnetic field and $\sigma_{T}$ is the Thomson cross
section. From the cooling time it is possible to calculate the cooling
coefficient, in fact $t_{\text{cool}}=p/\dot{p}$. Another effect influencing
the electron spectra is the spatial escape from the acceleration site:

$t_{\text{esc}}=\frac{R}{c}+\frac{R^{2}}{k_{\parallel}}$ (4)

where $k_{\parallel}=cr_{g}/9\zeta(k_{\text{res}})$ is the spatial diffusion
coefficient along the mean magnetic field, while $\zeta(k)=kW_{B}(k)/U_{B}$ is
the relative amplitude of the turbulent magnetic field energy density for a
given $k$. Note that if the mean free path is large (which implies a small
diffusion coefficient), the escape time is simply equal to the geometric
escape time, that is the ratio between the dimension of the acceleration
region and the particles velocity. Finally the last term of the equation (1)
describes the particle injection, identified here as the electrons accelerated
at the recollimation shock and advected downstream. Supposing a strong shock,
the injected electron density per unit time is equal to:

$I_{n}(\gamma)=I_{n,0}\,\gamma^{-2}\,e^{-\frac{\gamma}{\gamma_{\text{cut}}}}\quad\text{with}\quad\gamma_{\text{inj}}<\gamma<100\,\gamma_{\text{cut}}$
(5)

The extremes of the injection range are fixed and they are estimated thanks to
recent simulations of diffusive shock acceleration:
$\gamma_{\text{inj}}=10^{3}$ and $\gamma_{\text{cut}}=10^{5}$ (Zech & Lemoine,
2021). From the injected electron density per unit time it is possible to
obtain the injection distribution in the momentum space, in fact
$I_{n}(\gamma)=4\pi p^{2}m_{e}cI_{f}(p)$.

The second equation we have to consider describes in a phenomenological way
the time evolution of $W(k)$, the wavenumber energy distribution of the
turbulence, i.e. the energy density of waves with wavenumber between $k$ and
$k+dk$. Defining $Z(k)=W(k)/k^{2}$, in order to obtain a Fokker-Planck
equation in its canonical form, the equation can be written as (Eilek, 1979;
Zhou & Matthaeus, 1990; Kakuwa, 2016):

$\frac{\partial Z}{\partial t}=\frac{1}{k^{2}}\frac{\partial}{\partial
k}\left(k^{2}D_{k}\frac{\partial Z}{\partial
k}\right)+\frac{Z}{t_{\text{dam}}}+\frac{I_{W}}{k^{2}}$ (6)

where $D_{k}$ is the wavenumber diffusion coefficient, $t_{\text{dam}}$ the
damping time and $I_{W}$ the turbulence injection rate. We neglect spatial
effects, such as wave propagation outside the downstream region, and
dissipation mechanisms, such as Landau or cyclotron damping. Moreover there is
an important caveat: equation (6) describes an isotropic turbulent field,
while the acceleration time formula is calculated assuming that turbulence
consists exclusively in waves parallel or antiparallel to the magnetic field.
Even though it would be more correct to recompute equation (2) in the case of
isotropic turbulence, we adopt the convention used in the literature.

The wavenumber diffusion coefficient depends on which scenario we are
considering. Since we suppose that the jet magnetization is weak and the
turbulence is strong, the Kolmogorov phenomenology must be applied (Zhou &
Matthaeus, 1990; Miller & Roberts, 1995), which implies:

$t_{\text{cas}}=\frac{1}{k\beta_{w}c}\sqrt{\frac{2U_{B}}{kW(k)}}$ (7)

The damping time, instead, can be calculated from the acceleration term in
equation (1), requiring energy conservation. In fact the energy required to
accelerate and diffuse the electrons in the momentum space cannot be assumed
infinite, since in the considered scenario the jet magnetization is weak and
then the turbulence is extremely damped:

$t_{\text{dam}}=\left|\frac{4\pi
e^{2}\beta_{w}^{2}}{m_{e}ck}\int_{\gamma_{\text{res}}}\gamma^{2}\frac{\partial}{\partial\gamma}\left(\frac{n(\gamma)}{\gamma^{2}}\right)d\gamma\right|^{-1}$
(8)

where $\gamma_{\text{res}}$ is the resonant Lorentz factor, defined by
$k=2\pi/r_{g}(\gamma_{\text{res}})$. In fact an Alfvén wave with a fixed
wavenumber can interact only with electrons whose Larmor radius is greater
than its wavelength.

Finally the last term in equation (6) represents the turbulence injection. We
suppose that the dimension of the initial turbulence wavelength is in the
order of $L=R/10$:

$I_{W}(k)=I_{W,0}\,\delta(k-k_{0})$ (9)

where $\delta$ is the Dirac function and $k_{0}=2\pi/L$. The radiative
emission is calculated through the Synchrotron Self Compton model (see
Chiaberge & Ghisellini 1999 for the synchrotron emission and Blumenthal &
Gould 1970 for the inverse Compton emission), while for the beaming we use the
standard formula (see the discussion in Zech & Lemoine 2021).

##  3 Numerical method

For solving the two coupled non linear Fokker-Planck equations, we decided to
use the robust implicit Chang-Cooper algorithm (Chang & Cooper, 1970), in
which the diffusion and the cooling terms are regrouped in a single flux term,
whereas for the injection and the escape/damping terms we follow the technique
already implemented in other papers (Park & Petrosian, 1996; Chiaberge &
Ghisellini, 1999).

The time mesh is equally spaced and consists of $500$ points, while momentum
and wavenumber are divided in two logarithmic equally spaced meshes of $400$
points: the momentum mesh starts from the momentum corresponding to
$\gamma_{\text{min}}=10^{1}$ and ends when $\gamma_{\text{max}}=10^{9}$,
whereas the wavenumber space ranges from the injection wavenumber
($k_{\text{min}}=k_{0}$) to one thousand times the maximum resonant wavenumber
($k_{\text{max}}=1000\,k_{\text{res},\text{max}}$). Since the momentum and the
wavenumber meshes are logarithmic, we substitute both variables with the
corresponding logarithms, so both meshes result equally spaced with respect to
the new variables (Miller & Roberts, 1995). For the same reason, the
coefficient at midpoints are evaluated through the geometric mean, instead of
the usual arithmetic mean (Park & Petrosian, 1996).

The sizes of the Lorentz factor and wavenumber spaces are strongly linked to
the boundary conditions we imposed. The extremes of the Lorentz factor space
are chosen to be far from the range responsible for the radiative emission,
where the particle are the most, therefore we set
$f(p_{\text{min}})=f(p_{\text{max}})=0$. Considering the wavenumber space, we
set the flux through the lower limit to zero ($F_{W}(k_{\text{min}})=0$),
because the turbulence energy can not spread to wavenumber corresponding to
larger vortex, while the upper limit condition is less obvious, in fact both
setting the flux or the function itself to zero are unphysical conditions
(Park & Petrosian, 1996). For this reason we choose to extend the upper limit
of the wavenumber space beyond the maximum resonant wavenumber. In fact
between $k_{\text{res},\text{max}}$ and $k_{\text{max}}$ there are neither
source term nor damping term, therefore in this interval we have $W(k)\propto
k^{-5/3}$ and it is possible to obtain the value of the turbulence spectrum in
$k_{\text{max}}$ simply extending the function.

For the initial conditions, we suppose that at start there are no non-thermal
particles in the downstream region ($f(t_{0},p)=0$), while for the turbulence
spectrum we presume that the system is in equilibrium ($W(t_{0},k)\propto
k^{-5/3}$) and the normalization of the initial spectrum is supposed to be
much smaller than the injection one. A no null initial condition is necessary
for the injection to spread to higher wavenumbers, as it is clear from
equation (7).

Another problem for the system resolution is the non linearity of equation
(6), due to the presence of $W(k)$ in cascade time definition. Considering
equation (7), there are two possible solutions: the semi-implicit method, in
which we substitute $W(k)$ with the spectrum obtained in the previous time
step, or the fully implicit method, where $W(k)$ is set equal to the spectrum
we want to obtain in the considered time step. The first technique requires to
solve a linear system, but several time steps are necessary to converge, on
the other hand the second method converges despite of the width of the time
step, but it requires the resolution of a non linear system. We opt for the
second method, but we use an iterative procedure that requires solving several
linear systems: firstly we set $W(k)$ equal to the result obtained at the
previous time step, then we apply Chang Cooper algorithm, we substitute $W(k)$
with the solution, we reapply the algorithm and compare the two results. This
process continues until a defined convergence criterion is satisfied (Larsen
et al., 1985): in our case we require the relative error of $W(k)$ to be
smaller than $10^{-3}$.

##  4 Results

With the setup described above, our model depends on five free parameters: the
downstream region radius $R$, the Alfvén velocity $v_{a}$, the mean magnetic
field $B$, the injected powers of the electrons and the turbulence in the jet
frame, respectively $P_{n}^{\prime}$ and $P_{W}^{\prime}$. The latter two
quantities are defined as follows:

$\displaystyle P_{n}^{\prime}=V\int\gamma m_{e}c^{2}I_{n}(\gamma)\,d\gamma$
(10) $\displaystyle P_{W}^{\prime}=V\int I_{W}(k)\,dk$ (11)

where $V=10\pi R^{3}$ is the volume of the emitting region. From these
definitions we can calculate both the injection normalizations $I_{n,0}$ and
$I_{W,0}$.

We apply our model to the prototypical extreme TeV BL Lac object, i.e. 1ES
0229+200 ($z=0.1396$).

Figure 1: SED of 1ES 0229+200 (points with error bars, from Costamante et al.
2018) reproduced with the model presented in the paper (red solid line)

Fig. 1 shows a model realization ($R=1.2\times 10^{16}\,$\mathrm{cm}$$,
$v_{a}=2\times 10^{9}\,$\mathrm{c}\mathrm{m}\mathrm{/}\mathrm{s}$$,
$B=15.9\,$\mathrm{mG}$$, $P_{n}^{\prime}=7\times
10^{39}\,$\mathrm{erg}\mathrm{/}\mathrm{s}$$, $P_{W}^{\prime}=7\times
10^{39}\,$\mathrm{erg}\mathrm{/}\mathrm{s}$$ and with a relativistic Doppler
factor $\delta=25$, e.g. Tavecchio et al. 2022), which is in good agreement
with the data: in fact the X-ray peak is well reproduced and the model
realization presents a hard sub TeV spectrum with a gamma peak beyond
$1\,$\mathrm{TeV}$$. The free parameters are consistent with what we expect
from the theory. Since the emission is due to the post-reconfinement flow, $R$
is much smaller than the ordinary jet radius (Bodo & Tavecchio, 2018), whereas
the Alfvén velocity is large with respect to the usual values in order to make
the stochastic acceleration energy gain relevant for the non-thermal
particles. Moreover the magnetic field is in the order of the values obtained
by the phenomenological models that leave the acceleration process unspecified
(Biteau et al., 2020). Finally the injected power of the non-thermal particle
in the jet frame is in agreement with the theory predictions (Ghisellini et
al., 2010). To check the consistency of the model realization shown in Figure
1 further, we also calculate some derived quantities, i.e. the jet
magnetization $\sigma=4.4\times 10^{-3}$, the integrated relative amplitude of
the turbulent magnetic field $\zeta=0.1$, the background plasma number density
$n_{p}=2.8\ $\mathrm{c}\mathrm{m}^{-3}$$ and the final electron number density
$n_{e}=2.6\times 10^{-2}\ $\mathrm{c}\mathrm{m}^{-3}$$: they all assume
satisfactory values.

Figure 2: The time evolution of the electron energy distribution used to
produce Figure 1. The dashed lines represent different time steps (each equal
to $R/c$), while the solid line is the final distribution (reached at
$t=10R/c$). The time evolution is represented thanks to the red gradient. The
injection is not visible.

In Fig. 2 the time evolution of the electron energy distribution is
represented and several features can be noted. Firstly the maximum energy of
the distribution increases with time, reaching after $10\,R/c$ the equilibrium
value at $\gamma\sim 10^{6}$, where the acceleration and cooling time are
nearly equal (see Figure 4). Beyond the peak the final spectrum presents a
steep cut-off due to the synchrotron cooling, in fact in this range
$t_{\text{cool}}\ll t_{\text{acc}}$. Furthermore after some time the
distribution plunges for Lorentz factor smaller than $\gamma\sim 10^{4}$,
because for these particles the turbulence becomes negligible, than the escape
time is simply equal to $R/c$, as it is clear looking to the escape time in
Fig. 4.

Figure 3: The time evolution of the turbulence spectrum used to produce Figure
1. The dashed lines represent different time steps (each equal to $R/c$),
while the solid line is the final distribution. The time evolution is
represented thanks to the blue gradient.

Considering the time evolution of the turbulence spectrum in Fig. 3, initially
the injection cascades to larger wavenumbers, forming the usual power law
spectra, i.e. $W(k)\propto k^{-5/3}$. After spreading in all wavenumber space,
the non-thermal particles start to absorb the turbulence energy. For this
reason at wavenumbers larger than $k\sim
10^{-10}\,$\mathrm{c}\mathrm{m}^{-1}$$, i.e. when $t_{\text{dam}}\ll
t_{\text{cas}}$ (see Fig.4), the final turbulence spectrum is extremely
damped, while at shorter wavenumbers it remains a power law, in fact in this
range $t_{\text{cas}}\ll t_{\text{dam}}$.

Figure 4: Timescales characterizing the system at $t=10R/c$. Times related to
non-thermal particles are represented using solid lines, while for timescales
for turbulence we use dashed lines.

In order to better understand the role of the damping in shaping the electron
energy distribution and the corresponding produced SED, we decide to compare
our model with the data of 1ES0229+200, but excluding the damping term.

Figure 5: The time evolution of the electron energy distribution used to
produce Figure 6. The dashed lines represent different time steps (each equal
to $R/c$), while the solid line is the final distribution (reached at
$t=10R/c$). The time evolution is represented thanks to the red gradient. The
injection is not visible.

The absence of damping implies that a large amount of energy is stored in
waves at large wavenumbers ($k>10^{-10}\,$\mathrm{c}\mathrm{m}^{-1}$$, above
the cut-off resulting from damping visible in Fig. 3), causing the efficient
acceleration of particles at energies $\gamma\lesssim 10^{5}$. The rapid
acceleration, in turn, results in an electron energy distribution quite hard
before the peak (compare Figs. 2 and 5). Instead near and beyond the peak the
electron spectrum remains substantially unchanged: in fact even with damping,
in the corresponding wavenumber range the turbulence spectrum is the standard
power law.

The very narrow, maxwellian-like electron distribution translates into a very
narrow and hard synchrotron and SSC peaks in the SED. As shown in Fig. 6,
while the first peak is well reproduced, the rise of the SSC component in the
sub-TeV range is too hard with respect to the data, as already noted in
Tavecchio et al. (2022). Therefore, the inclusion of the damping, besides
making the model self-consistent, results in a softer gamma-ray spectrum, in
better agreement with the data.

Figure 6: SED of 1ES 0229+200 (points with error bars, from Costamante et al.
2018) reproduced with the model presented in the paper (red solid line)
excluding damping. Note the hard SSC spectrum.

##  5 Conclusions

In this paper we have presented a self-consistent model for extreme TeV BL Lac
objects, improving the simple version of Tavecchio et al. (2022). In
particular, the inclusion of damping for the turbulence results in an electron
distribution with an excess of particles below the peak, around $\gamma\simeq
10^{4}-10^{5}$, which produces a softer SSC component in the VHE band, in
better agreement with the data with respect to the simple model without wave
damping.

Our model is able to reproduce the SED of the prototypical extreme TeV BL Lac
object 1ES 0229+200, but there are several caveats. First of all it will be
necessary to compare our model with the data of other extreme TeV BL Lacs in
order to check if this hybrid emission scenario can be generalized to the
entire population. Since the classification of this type of sources is still
unclear, a different emission scenario could trigger a new way to categorize
these objects. Moreover the parameters obtained by the comparison between SEDs
and the model have to be checked with the results of MHD simulation. The
papers already published shows that the low magnetization scenario favors the
development of turbulence after the recollimation shock instead of the
formation of multiple shocks. A better characterization of the emission region
and turbulence based on dedicated MHD simulations will be presented in the
near future (Costa et al., in prep).

Considering the model in this paper, there are some possible improvements.
Firstly the Inverse Compton (i.e. SSC) cooling term could be added in the
diffusion equation of non-thermal electrons, which would make the equation (1)
non-linear too: if necessary, the same solution proposed for the non-linearity
of the diffusion equation of turbulence could be used. Moreover the
computational scheme used to solve the system in this paper could be replaced
with more complex and accurate algorithms, see for example Kundu et al. 2021
for an approach based on the employment of Runge-Kutta Implicit-Explicit
schemes.

In the following years the number of extreme TeV BL Lacs will probably
increase and their SED measurements will improve. In fact many data are going
to be collected, since these sources are already studied by several telescopes
(Costamante et al., 2018; Acciari et al., 2020) and they will also be objects
of interest for future gamma-ray facilities, such as CTA, because of their
extreme high energy emission. Finally, another observational feature useful to
test our model will be the polarization properties of the extreme TeV BL Lacs.
Polarimetric measures in the optical band are extremely difficult for this
type of sources because of the weak emission of the jet in this band,
dominated by the emission of the host galaxy. Instead in the X-ray band the
emission is mainly due to the synchrotron cooling of the jet particles
emitting in a highly turbulent field, therefore we expect a low degree of
polarization. This prediction could be proved by the IXPE mission (Weisskopf
et al., 2016), the first satellite able to measure the polarization in the
X-ray band, launched at the end of 2021.

## Data availability

Data available on request.

## References

  * Acciari et al. (2020) Acciari V. A., et al., 2020, ApJS, 247, 16
  * Aharonian et al. (2008) Aharonian F. A., Khangulyan D., Costamante L., 2008, MNRAS, 387, 1206
  * Biteau et al. (2020) Biteau J., et al., 2020, Nature Astronomy, 4, 124
  * Blandford & Eichler (1987) Blandford R., Eichler D., 1987, Phys. Rep., 154, 1
  * Blandford et al. (2019) Blandford R., Meier D., Readhead A., 2019, ARA&A, 57, 467
  * Blumenthal & Gould (1970) Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics, 42, 237
  * Bodo & Tavecchio (2018) Bodo G., Tavecchio F., 2018, A&A, 609, A122
  * Böttcher et al. (2008) Böttcher M., Dermer C. D., Finke J. D., 2008, ApJ, 679, L9
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Bresci et al. (2022) Bresci V., Lemoine M., Gremillet L., Comisso L., Sironi L., Demidem C., 2022, arXiv e-prints, p. arXiv:2206.08380
  * Chang & Cooper (1970) Chang J. S., Cooper G., 1970, Journal of Computational Physics, 6, 1
  * Chiaberge & Ghisellini (1999) Chiaberge M., Ghisellini G., 1999, MNRAS, 306, 551
  * Costamante et al. (2001) Costamante L., et al., 2001, A&A, 371, 512
  * Costamante et al. (2018) Costamante L., Bonnoli G., Tavecchio F., Ghisellini G., Tagliaferri G., Khangulyan D., 2018, MNRAS, 477, 4257
  * Eilek (1979) Eilek J. A., 1979, ApJ, 230, 373
  * Essey & Kusenko (2010) Essey W., Kusenko A., 2010, Astroparticle Physics, 33, 81
  * Fichet de Clairfontaine et al. (2021) Fichet de Clairfontaine G., Meliani Z., Zech A., Hervet O., 2021, A&A, 647, A77
  * Ghisellini et al. (2010) Ghisellini G., Tavecchio F., Foschini L., Ghirlanda G., Maraschi L., Celotti A., 2010, MNRAS, 402, 497
  * Ghisellini et al. (2017) Ghisellini G., Righi C., Costamante L., Tavecchio F., 2017, MNRAS, 469, 255
  * Gourgouliatos & Komissarov (2018) Gourgouliatos K. N., Komissarov S. S., 2018, Nature Astronomy, 2, 167
  * Kakuwa (2016) Kakuwa J., 2016, ApJ, 816, 24
  * Kakuwa et al. (2015) Kakuwa J., Toma K., Asano K., Kusunose M., Takahara F., 2015, MNRAS, 449, 551
  * Kundu et al. (2021) Kundu S., Vaidya B., Mignone A., 2021, ApJ, 921, 74
  * Larsen et al. (1985) Larsen E. W., Levermore C. D., Pomraning G. C., Sanderson J. G., 1985, Journal of Computational Physics, 61, 359
  * Lefa et al. (2011) Lefa E., Rieger F. M., Aharonian F., 2011, ApJ, 740, 64
  * Matsumoto et al. (2021) Matsumoto J., Komissarov S. S., Gourgouliatos K. N., 2021, MNRAS, 503, 4918
  * Miller & Roberts (1995) Miller J. A., Roberts D. A., 1995, ApJ, 452, 912
  * Park & Petrosian (1996) Park B. T., Petrosian V., 1996, ApJS, 103, 255
  * Romero et al. (2017) Romero G. E., Boettcher M., Markoff S., Tavecchio F., 2017, Space Sci. Rev., 207, 5
  * Sironi & Spitkovsky (2011) Sironi L., Spitkovsky A., 2011, ApJ, 726, 75
  * Sironi et al. (2015) Sironi L., Petropoulou M., Giannios D., 2015, MNRAS, 450, 183
  * Tavecchio et al. (2009) Tavecchio F., Ghisellini G., Ghirlanda G., Costamante L., Franceschini A., 2009, MNRAS, 399, L59
  * Tavecchio et al. (2022) Tavecchio F., Costa A., Sciaccaluga A., 2022, MNRAS, in press, p. arXiv:2207.12766
  * Weisskopf et al. (2016) Weisskopf M. C., et al., 2016, in den Herder J.-W. A., Takahashi T., Bautz M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray. p. 990517
  * Zech & Lemoine (2021) Zech A., Lemoine M., 2021, A&A, 654, A96
  * Zhou & Matthaeus (1990) Zhou Y., Matthaeus W. H., 1990, J. Geophys. Res., 95, 14881


parameters, rather than via optimistic value estimates, performs nearly as
well as ReLOAD. This is likely because the CMDP only has two states, so using
the gradient computed from value estimates obtained from an old batch of data
will likely be close to that obtained from the current batch of data, as
there’s more likely to be significant overlap. If we look to higher-
dimensional tasks like control suite (Fig. 7), there is a much larger gap
between the direct OGD approach and ReLOAD as the state-action pairs used for
value estimation are more likely to differ from one minibatch to another.

Figure 11: Learning curves for $\mu$-IMPALA and ReLOAD-IMPALA on Walker, Walk
with a height constraint (left) and Reacher, Easy with a velocity constraint
(right). The horizontal dotted line indicates the value of the constraint. In
each case, ReLOAD-IMPALA significantly dampens oscillations compared to
$\mu$-IMPALA.

Figure 12: Learning curves for $\mu$-DMPO and ReLOAD-DMPO on Walker, Walk with a velocity constraint (left), Quadruped, Walk with a torque constraint (center), and Humanoid, Walk with a height constraint (right). The horizontal dotted line indicates the value of the constraint. In each case, ReLOAD-DMPO significantly dampens oscillations compared to $\mu$-DMPO. Domain | Algorithm | Safety-Coeff/Threshold | Weighted Reward | $\mathbf{R_{penalized}}$ | Task Reward |  Constraint Violation
---|---|---|---|---|---|---
Walker | 0.1-D4PG | 0.05/0.057 | $-153.8\pm 12.1$ | $552.8\pm 12.0$ | $979.6\pm 0.3$ | $426.8\pm 12.1$
|  | 0.05/0.077 | $-175.4\pm 20.0$ | $544.3\pm 19.9$ | $979.1\pm 0.5$ | $434.8\pm 20.0$
|  | 0.05/0.097 | $-75.1\pm 20.5$ | $581.8\pm 20.3$ | $978.6\pm 0.4$ | $396.8\pm 20.5$
|  | 0.1/0.057 | $533.9\pm 9.0$ | $811.3\pm 8.7$ | $978.9\pm 0.4$ | $167.6\pm 9.0$
|  | 0.1/0.077 | $630.2\pm 7.9$ | $847.7\pm 7.7$ | $979.1\pm 0.3$ | $131.4\pm 7.9$
|  | 0.1/0.097 | $652.4\pm 5.8$ | $855.8\pm 5.7$ | $978.7\pm 0.3$ | $122.9\pm 5.8$
|  | 0.2/0.057 | $982.0\pm 0.2$ | $982.0\pm 0.2$ | $982.0\pm 0.2$ | $0.0\pm 0.0$
|  | 0.2/0.077 | $981.8\pm 0.2$ | $981.8\pm 0.2$ | $981.8\pm 0.2$ | $0.0\pm 0.0$
|  | 0.2/0.097 | $981.5\pm 0.3$ | $981.5\pm 0.3$ | $981.5\pm 0.3$ | $0.0\pm 0.0$
|  | 0.3/0.057 | $982.7\pm 0.3$ | $982.7\pm 0.3$ | $982.7\pm 0.3$ | $0.0\pm 0.0$
|  | 0.3/0.077 | $981.8\pm 0.3$ | $981.8\pm 0.3$ | $981.8\pm 0.3$ | $0.0\pm 0.0$
|  | 0.3/0.097 | $982.0\pm 0.1$ | $982.0\pm 0.1$ | $982.0\pm 0.1$ | $0.0\pm 0.0$
| RC-D4PG | 0.05/0.057 | $268.7\pm 34.6$ | $296.1\pm 28.3$ | $312.7\pm 34.0$ | $16.6\pm 6.4$
|  | 0.05/0.077 | $325.9\pm 44.4$ | $346.6\pm 36.9$ | $359.0\pm 43.8$ | $12.5\pm 7.2$
|  | 0.05/0.097 | $329.1\pm 33.4$ | $334.8\pm 32.8$ | $338.2\pm 33.4$ | $3.4\pm 1.3$
|  | 0.1/0.057 | $914.3\pm 10.9$ | $915.0\pm 10.8$ | $915.5\pm 10.9$ | $0.5\pm 0.4$
|  | 0.1/0.077 | $942.7\pm 1.6$ | $945.4\pm 1.2$ | $947.1\pm 1.4$ | $1.7\pm 0.8$
|  | 0.1/0.097 | $956.1\pm 3.3$ | $959.3\pm 3.3$ | $961.3\pm 3.2$ | $2.0\pm 0.8$
|  | 0.2/0.057 | $980.5\pm 0.3$ | $981.0\pm 0.4$ | $981.4\pm 0.2$ | $0.3\pm 0.2$
|  | 0.2/0.077 | $979.1\pm 1.0$ | $980.7\pm 1.1$ | $981.7\pm 0.3$ | $1.0\pm 1.0$
|  | 0.2/0.097 | $977.7\pm 0.8$ | $980.3\pm 0.9$ | $981.9\pm 0.3$ | $1.6\pm 0.8$
|  | 0.3/0.057 | $978.4\pm 0.5$ | $980.5\pm 0.5$ | $981.8\pm 0.3$ | $1.3\pm 0.4$
|  | 0.3/0.077 | $978.6\pm 0.9$ | $980.5\pm 0.7$ | $981.6\pm 0.3$ | $1.1\pm 0.8$
|  | 0.3/0.097 | $979.6\pm 0.6$ | $980.9\pm 0.6$ | $981.7\pm 0.4$ | $0.8\pm 0.5$
| MetaL | 0.05/0.057 | $-168.5\pm 10.5$ | $540.0\pm 9.6$ | $968.0\pm 1.5$ | $428.0\pm 10.4$
|  | 0.05/0.077 | $-138.8\pm 13.7$ | $551.7\pm 13.5$ | $968.9\pm 1.4$ | $417.1\pm 13.6$
|  | 0.05/0.097 | $-131.1\pm 13.0$ | $553.9\pm 12.2$ | $967.8\pm 1.0$ | $413.8\pm 13.0$
|  | 0.1/0.057 | $663.4\pm 19.8$ | $853.7\pm 19.8$ | $968.6\pm 1.2$ | $114.9\pm 19.8$
|  | 0.1/0.077 | $627.4\pm 17.8$ | $841.7\pm 17.1$ | $971.2\pm 1.1$ | $129.5\pm 17.8$
|  | 0.1/0.097 | $766.9\pm 19.8$ | $893.3\pm 19.7$ | $969.7\pm 0.7$ | $76.4\pm 19.7$
|  | 0.2/0.057 | $981.3\pm 0.5$ | $982.0\pm 0.6$ | $982.5\pm 0.2$ | $0.4\pm 0.4$
|  | 0.2/0.077 | $980.2\pm 0.3$ | $981.5\pm 0.3$ | $982.3\pm 0.1$ | $0.8\pm 0.3$
|  | 0.2/0.097 | $978.6\pm 0.8$ | $980.7\pm 0.7$ | $981.9\pm 0.3$ | $1.3\pm 0.7$
|  | 0.3/0.057 | $977.7\pm 0.9$ | $980.6\pm 1.0$ | $982.3\pm 0.4$ | $1.7\pm 0.8$
|  | 0.3/0.077 | $978.4\pm 0.5$ | $980.6\pm 0.6$ | $981.9\pm 0.2$ | $1.3\pm 0.5$
|  | 0.3/0.097 | $979.6\pm 0.5$ | $980.9\pm 0.5$ | $981.6\pm 0.2$ | $0.8\pm 0.5$
| ReLOAD | 0.05/0.057 | $367.0\pm 114.8$ | $626.9\pm 113.4$ | $783.8\pm 96.9$ | $157.0\pm 9.6$
|  | 0.05/0.077 | $423.3\pm 80.1$ | $716.7\pm 43.5$ | $893.9\pm 26.1$ | $177.2\pm 24.0$
|  | 0.05/0.097 | $315.8\pm 43.1$ | $677.7\pm 28.9$ | $896.4\pm 22.1$ | $218.6\pm 10.3$
|  | 0.1/0.057 | $754.1\pm 39.0$ | $875.4\pm 23.7$ | $948.7\pm 15.8$ | $73.3\pm 10.1$
|  | 0.1/0.077 | $717.9\pm 26.8$ | $877.1\pm 11.2$ | $973.3\pm 4.8$ | $96.2\pm 9.9$
|  | 0.1/0.097 | $746.9\pm 28.4$ | $889.3\pm 12.0$ | $975.3\pm 4.7$ | $86.0\pm 10.4$
|  | 0.2/0.057 | $907.7\pm 6.9$ | $950.6\pm 5.0$ | $976.6\pm 4.1$ | $26.0\pm 1.9$
|  | 0.2/0.077 | $881.5\pm 9.7$ | $940.7\pm 5.9$ | $976.4\pm 4.3$ | $35.7\pm 3.1$
|  | 0.2/0.097 | $882.7\pm 7.8$ | $942.1\pm 5.2$ | $977.9\pm 4.0$ | $35.9\pm 2.4$
|  | 0.3/0.057 | $938.3\pm 5.8$ | $963.0\pm 5.0$ | $977.9\pm 4.3$ | $14.9\pm 1.2$
|  | 0.3/0.077 | $930.4\pm 5.6$ | $960.5\pm 4.6$ | $978.8\pm 4.0$ | $18.2\pm 1.3$
|  | 0.3/0.097 | $925.8\pm 5.8$ | $958.3\pm 4.7$ | $977.9\pm 4.0$ | $19.6\pm 1.5$
Table 2: Mean performance for RWRL-Walker averaged over 8 random seeds for each of the 12 constraint settings in the RWRL suite, where lower values of the safety coefficient and threshold indicate more challenging tasks. $\pm$ values denote one standard error.  Domain | Algorithm | Safety-Coeff/Threshold | Weighted Reward | $\mathbf{R_{penalized}}$ | Task Reward |  Constraint Violation
---|---|---|---|---|---|---
Cartpole | 0.1-D4PG | 0.05/0.07 | $833.8\pm 1.3$ | $352.6\pm 1.2$ | $851.7\pm 0.2$ | $499.1\pm 1.3$
|  | 0.05/0.09 | $834.3\pm 1.5$ | $372.6\pm 1.5$ | $851.5\pm 0.1$ | $478.9\pm 1.5$
|  | 0.05/0.115 | $835.0\pm 0.9$ | $399.7\pm 0.8$ | $851.2\pm 0.1$ | $451.5\pm 0.9$
|  | 0.1/0.07 | $844.9\pm 0.7$ | $568.2\pm 0.7$ | $855.2\pm 0.1$ | $287.0\pm 0.7$
|  | 0.1/0.09 | $845.5\pm 1.0$ | $588.3\pm 1.0$ | $855.1\pm 0.2$ | $266.8\pm 1.0$
|  | 0.1/0.115 | $846.4\pm 0.9$ | $612.9\pm 0.9$ | $855.1\pm 0.2$ | $242.2\pm 0.9$
|  | 0.2/0.07 | $850.7\pm 4.8$ | $763.1\pm 6.3$ | $853.9\pm 2.4$ | $90.8\pm 4.2$
|  | 0.2/0.09 | $854.7\pm 1.3$ | $791.3\pm 1.5$ | $857.0\pm 0.5$ | $65.7\pm 1.2$
|  | 0.2/0.115 | $852.0\pm 4.4$ | $807.6\pm 5.8$ | $853.6\pm 2.4$ | $46.0\pm 3.7$
|  | 0.3/0.07 | $857.5\pm 0.3$ | $825.5\pm 0.3$ | $858.7\pm 0.1$ | $33.3\pm 0.3$
|  | 0.3/0.09 | $858.4\pm 0.6$ | $845.4\pm 0.5$ | $858.9\pm 0.1$ | $13.5\pm 0.5$
|  | 0.3/0.115 | $856.7\pm 1.9$ | $856.6\pm 2.0$ | $856.7\pm 1.9$ | $0.1\pm 0.1$
| RC-D4PG | 0.05/0.07 | $278.6\pm 114.7$ | $158.7\pm 21.7$ | $283.1\pm 88.2$ | $124.4\pm 73.3$
|  | 0.05/0.09 | $260.4\pm 84.4$ | $144.8\pm 43.9$ | $264.7\pm 66.4$ | $119.9\pm 52.2$
|  | 0.05/0.115 | $422.8\pm 105.5$ | $269.1\pm 23.0$ | $428.5\pm 83.5$ | $159.4\pm 64.4$
|  | 0.1/0.07 | $262.6\pm 48.7$ | $180.4\pm 49.6$ | $265.7\pm 43.1$ | $85.3\pm 22.6$
|  | 0.1/0.09 | $534.4\pm 104.0$ | $296.2\pm 93.7$ | $543.3\pm 81.9$ | $247.1\pm 64.1$
|  | 0.1/0.115 | $710.0\pm 45.1$ | $537.5\pm 35.5$ | $716.4\pm 38.8$ | $178.9\pm 23.0$
|  | 0.2/0.07 | $310.9\pm 90.6$ | $298.3\pm 86.4$ | $311.4\pm 90.3$ | $13.1\pm 7.5$
|  | 0.2/0.09 | $770.0\pm 18.0$ | $760.6\pm 16.4$ | $770.3\pm 17.7$ | $9.7\pm 3.0$
|  | 0.2/0.115 | $847.7\pm 2.4$ | $847.4\pm 2.4$ | $847.7\pm 2.4$ | $0.3\pm 0.1$
|  | 0.3/0.07 | $845.2\pm 3.4$ | $844.7\pm 3.6$ | $845.3\pm 3.4$ | $0.6\pm 0.3$
|  | 0.3/0.09 | $857.3\pm 0.2$ | $857.1\pm 0.1$ | $857.3\pm 0.1$ | $0.2\pm 0.1$
|  | 0.3/0.115 | $859.2\pm 0.3$ | $858.7\pm 0.2$ | $859.2\pm 0.1$ | $0.4\pm 0.3$
| MetaL | 0.05/0.07 | $832.9\pm 1.1$ | $348.2\pm 1.1$ | $850.9\pm 0.2$ | $502.7\pm 1.1$
|  | 0.05/0.09 | $829.5\pm 2.3$ | $365.1\pm 2.5$ | $846.8\pm 2.1$ | $481.6\pm 1.0$
|  | 0.05/0.115 | $834.4\pm 1.7$ | $391.9\pm 1.6$ | $850.9\pm 0.2$ | $458.9\pm 1.7$
|  | 0.1/0.07 | $843.4\pm 1.0$ | $580.3\pm 0.9$ | $853.2\pm 0.2$ | $272.9\pm 1.0$
|  | 0.1/0.09 | $844.3\pm 2.3$ | $599.3\pm 2.3$ | $853.4\pm 0.2$ | $254.1\pm 2.3$
|  | 0.1/0.115 | $845.5\pm 1.2$ | $623.7\pm 1.3$ | $853.7\pm 0.1$ | $230.0\pm 1.2$
|  | 0.2/0.07 | $851.1\pm 2.6$ | $787.3\pm 3.3$ | $853.5\pm 2.0$ | $66.2\pm 1.6$
|  | 0.2/0.09 | $850.8\pm 3.8$ | $809.3\pm 4.2$ | $852.3\pm 2.1$ | $43.1\pm 3.1$
|  | 0.2/0.115 | $855.1\pm 2.1$ | $831.8\pm 1.9$ | $855.9\pm 0.2$ | $24.1\pm 2.1$
|  | 0.3/0.07 | $856.1\pm 1.9$ | $839.9\pm 1.5$ | $856.7\pm 0.4$ | $16.8\pm 1.9$
|  | 0.3/0.09 | $852.8\pm 3.2$ | $848.5\pm 4.1$ | $852.9\pm 2.7$ | $4.4\pm 1.6$
|  | 0.3/0.115 | $857.3\pm 2.0$ | $856.9\pm 2.1$ | $857.3\pm 2.0$ | $0.4\pm 0.2$
| ReLOAD | 0.05/0.07 | $875.8\pm 2.6$ | $798.3\pm 9.4$ | $878.7\pm 2.4$ | $80.4\pm 7.4$
|  | 0.05/0.09 | $871.5\pm 3.0$ | $781.0\pm 11.2$ | $874.8\pm 2.8$ | $93.8\pm 8.7$
|  | 0.05/0.115 | $869.1\pm 3.0$ | $769.0\pm 11.3$ | $872.8\pm 2.7$ | $103.9\pm 8.9$
|  | 0.1/0.07 | $875.8\pm 2.6$ | $801.9\pm 9.0$ | $878.5\pm 2.4$ | $76.6\pm 6.9$
|  | 0.1/0.09 | $869.5\pm 2.9$ | $780.7\pm 9.9$ | $872.8\pm 2.7$ | $92.1\pm 7.6$
|  | 0.1/0.115 | $871.7\pm 3.0$ | $786.6\pm 10.1$ | $874.8\pm 2.8$ | $88.2\pm 7.7$
|  | 0.2/0.07 | $869.9\pm 2.9$ | $795.1\pm 8.1$ | $872.7\pm 2.7$ | $77.5\pm 6.0$
|  | 0.2/0.09 | $871.9\pm 3.0$ | $798.7\pm 8.9$ | $874.7\pm 2.8$ | $76.0\pm 6.5$
|  | 0.2/0.115 | $869.8\pm 2.9$ | $789.2\pm 8.6$ | $872.7\pm 2.7$ | $83.5\pm 6.4$
|  | 0.3/0.07 | $870.1\pm 2.8$ | $802.9\pm 7.5$ | $872.6\pm 2.7$ | $69.7\pm 5.5$
|  | 0.3/0.09 | $867.8\pm 2.6$ | $797.1\pm 7.3$ | $870.5\pm 2.5$ | $73.4\pm 5.5$
|  | 0.3/0.115 | $867.6\pm 2.6$ | $793.3\pm 7.3$ | $870.4\pm 2.4$ | $77.1\pm 5.5$
Table 3: Mean performance for RWRL-Cartpole averaged over 8 random seeds for each of the 12 constraint settings in the RWRL suite, where lower values of the safety coefficient and threshold indicate more challenging tasks. $\pm$ values denote one standard error.  Domain | Algorithm | Safety-Coeff/Threshold | Weighted Reward | $\mathbf{R_{penalized}}$ | Task Reward |  Constraint Violation
---|---|---|---|---|---|---
Quadruped | 0.1-D4PG | 0.05/0.545 | $109.5\pm 0.0$ | $546.6\pm 0.0$ | $999.6\pm 0.0$ | $453.0\pm 0.0$
|  | 0.05/0.645 | $305.9\pm 0.0$ | $646.5\pm 0.0$ | $999.5\pm 0.0$ | $353.0\pm 0.0$
|  | 0.05/0.745 | $502.4\pm 0.0$ | $746.5\pm 0.0$ | $999.5\pm 0.0$ | $253.0\pm 0.0$
|  | 0.1/0.545 | $523.7\pm 0.0$ | $547.5\pm 0.0$ | $999.5\pm 0.0$ | $452.0\pm 0.0$
|  | 0.1/0.645 | $629.0\pm 0.0$ | $647.5\pm 0.0$ | $999.5\pm 0.0$ | $352.0\pm 0.0$
|  | 0.1/0.745 | $734.2\pm 0.0$ | $747.5\pm 0.0$ | $999.5\pm 0.0$ | $252.0\pm 0.0$
|  | 0.2/0.545 | $988.7\pm 56.4$ | $942.9\pm 56.3$ | $999.3\pm 0.0$ | $56.4\pm 56.4$
|  | 0.2/0.645 | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $0.0\pm 0.0$
|  | 0.2/0.745 | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $0.0\pm 0.0$
|  | 0.3/0.545 | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $0.0\pm 0.0$
|  | 0.3/0.645 | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $999.3\pm 0.0$ | $0.0\pm 0.0$
|  | 0.3/0.745 | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $0.0\pm 0.0$
| RC-D4PG | 0.05/0.545 | $-157.5\pm 71.1$ | $279.5\pm 71.1$ | $732.5\pm 71.1$ | $453.0\pm 0.0$
|  | 0.05/0.645 | $28.1\pm 60.0$ | $368.4\pm 60.0$ | $721.2\pm 60.0$ | $352.8\pm 0.2$
|  | 0.05/0.745 | $298.6\pm 49.0$ | $542.7\pm 49.0$ | $795.7\pm 49.0$ | $253.0\pm 0.0$
|  | 0.1/0.545 | $454.9\pm 116.0$ | $469.8\pm 158.4$ | $751.8\pm 81.5$ | $282.1\pm 82.6$
|  | 0.1/0.645 | $948.0\pm 43.8$ | $950.4\pm 44.7$ | $995.7\pm 1.2$ | $45.3\pm 43.8$
|  | 0.1/0.745 | $676.6\pm 56.6$ | $686.6\pm 73.0$ | $875.6\pm 38.7$ | $189.0\pm 41.2$
|  | 0.2/0.545 | $999.1\pm 0.0$ | $999.1\pm 0.0$ | $999.1\pm 0.0$ | $0.0\pm 0.0$
|  | 0.2/0.645 | $998.7\pm 1.4$ | $996.5\pm 1.4$ | $999.1\pm 0.1$ | $2.6\pm 1.4$
|  | 0.2/0.745 | $997.9\pm 4.1$ | $992.8\pm 4.1$ | $999.1\pm 0.1$ | $6.3\pm 4.1$
|  | 0.3/0.545 | $998.8\pm 2.6$ | $996.6\pm 2.6$ | $999.2\pm 0.0$ | $2.6\pm 2.6$
|  | 0.3/0.645 | $999.2\pm 0.2$ | $999.0\pm 0.2$ | $999.2\pm 0.0$ | $0.2\pm 0.2$
|  | 0.3/0.745 | $999.1\pm 0.0$ | $999.1\pm 0.0$ | $999.1\pm 0.0$ | $0.0\pm 0.0$
| MetaL | 0.05/0.545 | $108.6\pm 0.1$ | $545.6\pm 0.1$ | $998.6\pm 0.1$ | $452.9\pm 0.1$
|  | 0.05/0.645 | $305.3\pm 0.1$ | $645.8\pm 0.1$ | $998.7\pm 0.0$ | $352.9\pm 0.1$
|  | 0.05/0.745 | $501.6\pm 0.0$ | $745.7\pm 0.0$ | $998.7\pm 0.0$ | $253.0\pm 0.0$
|  | 0.1/0.545 | $523.0\pm 0.3$ | $546.8\pm 0.3$ | $998.4\pm 0.1$ | $451.6\pm 0.3$
|  | 0.1/0.645 | $674.4\pm 44.0$ | $690.7\pm 44.1$ | $998.7\pm 0.1$ | $308.0\pm 44.0$
|  | 0.1/0.745 | $766.7\pm 31.4$ | $778.3\pm 31.5$ | $998.4\pm 0.1$ | $220.1\pm 31.4$
|  | 0.2/0.545 | $998.3\pm 4.6$ | $994.2\pm 4.6$ | $999.2\pm 0.0$ | $4.9\pm 4.6$
|  | 0.2/0.645 | $999.0\pm 0.7$ | $998.4\pm 0.7$ | $999.1\pm 0.0$ | $0.7\pm 0.7$
|  | 0.2/0.745 | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $999.2\pm 0.0$ | $0.0\pm 0.0$
|  | 0.3/0.545 | $999.2\pm 0.3$ | $998.9\pm 0.3$ | $999.2\pm 0.0$ | $0.3\pm 0.3$
|  | 0.3/0.645 | $999.2\pm 0.2$ | $999.0\pm 0.2$ | $999.2\pm 0.0$ | $0.2\pm 0.2$
|  | 0.3/0.745 | $998.7\pm 1.8$ | $995.7\pm 1.8$ | $999.2\pm 0.0$ | $3.5\pm 1.8$
| ReLOAD | 0.05/0.545 | $651.8\pm 27.6$ | $657.7\pm 16.6$ | $663.8\pm 4.7$ | $6.1\pm 25.4$
|  | 0.05/0.645 | $657.7\pm 3.4$ | $659.3\pm 3.4$ | $661.0\pm 0.9$ | $1.7\pm 3.2$
|  | 0.05/0.745 | $650.0\pm 19.2$ | $654.7\pm 18.7$ | $659.5\pm 1.9$ | $4.9\pm 19.1$
|  | 0.1/0.545 | $996.4\pm 1.0$ | $996.4\pm 1.0$ | $997.1\pm 0.3$ | $0.7\pm 0.9$
|  | 0.1/0.645 | $993.4\pm 3.7$ | $993.7\pm 3.5$ | $998.4\pm 0.5$ | $4.7\pm 3.6$
|  | 0.1/0.745 | $982.3\pm 7.3$ | $983.2\pm 7.2$ | $999.1\pm 0.1$ | $16.0\pm 7.3$
|  | 0.2/0.545 | $998.0\pm 2.6$ | $992.6\pm 2.6$ | $999.3\pm 0.0$ | $6.7\pm 2.6$
|  | 0.2/0.645 | $997.4\pm 3.3$ | $988.9\pm 3.3$ | $999.3\pm 0.0$ | $10.4\pm 3.3$
|  | 0.2/0.745 | $996.2\pm 4.6$ | $982.8\pm 4.6$ | $999.3\pm 0.0$ | $16.5\pm 4.6$
|  | 0.3/0.545 | $998.4\pm 2.6$ | $992.5\pm 2.6$ | $999.3\pm 0.0$ | $6.8\pm 2.6$
|  | 0.3/0.645 | $998.0\pm 3.1$ | $989.7\pm 3.2$ | $999.3\pm 0.1$ | $9.6\pm 3.1$
|  | 0.3/0.745 | $997.9\pm 3.7$ | $988.8\pm 3.7$ | $999.3\pm 0.0$ | $10.6\pm 3.7$
Table 4: Mean performance for RWRL-Quadruped averaged over 8 random seeds for each of the 12 constraint settings in the RWRL suite, where lower values of the safety coefficient and threshold indicate more challenging tasks. $\pm$ values denote one standard error.  Domain (Constraint) | Algorithm | Weighted Reward | Task Reward |  Constraint Violation
---|---|---|---|---
Walker (Height) |  $\mu$-IMPALA | $354.5\pm 27.7$ | $360.6\pm 32.4$ | $5.5\pm 4.6$
| OGD-IMPALA | $487.5\pm 2.8$ | $510.1\pm 15.1$ | $20.4\pm 8.3$
|  $\mu^{\star}$-IMPALA | $356.7\pm 32.8$ | $527.5\pm 65.0$ | $153.6\pm 29.0$
| PID-IMPALA | $408.1\pm 12.7$ | $462.6\pm 27.4$ | $49.0\pm 19.8$
| RNTR-IMPALA | $386.6\pm 22.7$ | $399.0\pm 28.2$ | $11.1\pm 6.5$
| ReLOAD-IMPALA | $549.0\pm 2.0$ | $592.2\pm 16.4$ | $38.8\pm 10.8$
Reacher (Velocity) |  $\mu$-IMPALA | $436.9\pm 31.7$ | $450.1\pm 68.4$ | $117.3\pm 26.0$
| OGD-IMPALA | $22.4\pm 10.9$ | $44.6\pm 13.1$ | $197.7\pm 22.5$
|  $\mu^{\star}$-IMPALA | $893.6\pm 8.5$ | $913.6\pm 18.0$ | $178.1\pm 11.6$
| PID-IMPALA | $428.5\pm 25.6$ | $439.5\pm 72.0$ | $97.6\pm 22.2$
| RNTR-IMPALA | $841.7\pm 102.8$ | $850.0\pm 103.0$ | $73.7\pm 11.6$
| ReLOAD-IMPALA | $938.1\pm 2.9$ | $961.8\pm 11.7$ | $209.7\pm 14.1$
Table 5: Oscillating control suite results for IMPALA-based agents, averaged over 8 random seeds. $\pm$ values indicate one standard error. Domain (Constraint) | Algorithm | Weighted Reward | Task Reward |  Constraint Violation
---|---|---|---|---
Walker (Velocity) |  $\mu$-DMPO | $195.9\pm 57.6$ | $441.0\pm 57.5$ | $150.9\pm 2.1$
| OGD-DMPO | $108.2\pm 10.0$ | $357.5\pm 10.0$ | $153.5\pm 1.1$
|  $\mu^{\star}$-DMPO | $118.1\pm 15.9$ | $932.7\pm 3.8$ | $501.5\pm 15.4$
| PID-DMPO | $300.1\pm 18.6$ | $544.2\pm 18.5$ | $150.3\pm 1.6$
| ReLOAD-DMPO | $321.0\pm 49.3$ | $569.4\pm 49.2$ | $152.9\pm 2.0$
Quadruped (Torque) |  $\mu$-DMPO | $-254.8\pm 23.5$ | $567.1\pm 22.1$ | $558.3\pm 8.0$
| OGD-DMPO | $-220.1\pm 20.8$ | $664.6\pm 20.7$ | $601.0\pm 1.3$
|  $\mu^{\star}$-DMPO | $-328.4\pm 36.1$ | $200.4\pm 33.7$ | $359.2\pm 13.0$
| PID-DMPO | $-232.2\pm 21.1$ | $655.4\pm 20.7$ | $602.9\pm 4.2$
| ReLOAD-DMPO | $-137.0\pm 55.7$ | $768.4\pm 51.7$ | $615.0\pm 20.6$
Humanoid (Height) |  $\mu$-DMPO | $92.5\pm 76.8$ | $685.8\pm 59.4$ | $576.7\pm 48.7$
| OGD-DMPO | $1.9\pm 64.3$ | $606.1\pm 49.6$ | $587.2\pm 40.9$
|  $\mu^{\star}$-DMPO | $-45.9\pm 0.1$ | $0.0\pm 0.0$ | $44.7\pm 0.1$
| PID-DMPO | $39.9\pm 67.0$ | $538.3\pm 50.0$ | $484.4\pm 44.6$
| ReLOAD-DMPO | $114.4\pm 44.3$ | $650.1\pm 36.3$ | $520.6\pm 25.2$
Table 6: Oscillating control suite results for DMPO-based agents, averaged over 8 random seeds. $\pm$ values indicate one standard error. Hyperparameter |  Value
---|---
regularization weight $\alpha_{\Omega}$ | 1e-2
value loss weight $\alpha_{V}$ | 0.25
policy loss weight $\alpha_{\pi}$ | 1.0
discount factor $\gamma$ | 0.99
RMSProp decay | 0.99
RMSProp $\epsilon$ | 1e-4
initial step size | 6e-4
final step size | 6e-6
max gradient norm | 0.2
inner loop steps | 5
initial Bregman step size | 2.0
final Bregman step size | 0.5
initial $\mu$ learning rate | 1e-1
final $\mu$ learning rate | 1e-3
network hidden units per layer | 256
network depth | 3
training duration (environment steps) | 10e6
Table 7: Hyperparameter settings for IMPALA experiments. Hyperparameter |  Value
---|---
batch size | 256
replay size | 2e6
optimizer | Adam
regularization weight | 1e-2
value loss weight | 0.25
discount factor $\gamma$ | 0.99
initial step size | 3e-4
final step size | 1e-5
initial $\mu$ learning rate | 1e-1
final $\mu$ learning rate | 1e-3
network hidden units per layer | 256
network depth | 3
training duration | 15e6
Table 8: Hyperparameter settings for DMPO experiments.

Figure 13: ReLOAD-DMPO outperforms $\mu$-DMPO and matches the performance of
MO-DMPO on the RWRL Suite. Here, we trained agents on the comparatively easy
RWRL settings used by (Huang et al., 2022) over 8 random seeds. Error bars
denote one standard error.